Chronic infections reduce mortality following unrelated secondary infections in Drosophila melanogaster
Abigail Wukitch, Moria Chambers and Owais Gilani
6/17/2022
Q1 Does the protection conferred by primary infection with S. marcescens or E. faecalis depend on the infectious dose of the secondary infection with P. rettgeri?
Study design: 
- Previous work: It has been previously published that chronic infection with S. marcescens and E. faecalis provide significant protection against secondary infection with P. rettgeri initiated by an infection dose of 23nL of 0.1 Abs600nm (~8000 CFU/fly) (Chambers et al. 2019)
- We used infectious doses at least 10 fold above and below this dose, to assess how sensitive this protection is to dose.
- Flies were first infected with S. marcescens or E. faecalis and then a week later were given P. rettgeri infections with a range of doses (0.001-5.0 Abs600nm) and survival monitored for at least three days.
- Initially analyzed using Cox Proportional Hazards mixed effects model with date of injection, primary injection and secondary dose as fixed effects and random effect terms accounting for experimental structure (e.g. date of injection, housing of groups of flies in vials). However, this model violated the assumptions of Cox Proportional Hazards.
- Therefore we switched to using log-rank pairwise comparisons to
determine whether chronic infection significantly impacted survival of
secondary infection at an individual dose.
- LogRank comparisons were done on data from each date individually
- If experiments on all dates showed higher survival in conditions with chronic infection, p-values were combined using Fisher’s test and only the final p-value is reported in the text.
- If there was a variable effect of chronic infection, the number of individual experiments with significant effects is reported.
- Conclusions
- S. marcescens chronic infection provided significant protection during secondary infection with P. rettgeri at all infectious doses (p<0.001)
- E. faecalis chronic infection with offered significant protection against three out of the four doses (P < 0.05). However, the one dose that did not show significant protection (Abs600nm = 0.1) had a borderline p-value = 0.06
Figure 1 - S. marcescens Chronic Infection
- Mortality over 7 days following secondary infection with P. rettgeri in flies chronically infected with S. marcescens versus control flies not given a chronic infection.
Total Graph
#Sm_Pr$primary <- ordered(Sm_Pr$primary, levels = c("PBS", "S.m"))
table(Sm_Pr$primary, Sm_Pr$secondary_dose)##
## 0 0.01 0.1 1 2 3 5
## PBS 115 118 206 283 263 80 165
## S.m 114 119 196 275 238 81 159#Sm_bacterial_load$primary <- ordered(Sm_bacterial_load$primary, levels = c("PBS", "S.m"))
#table(Sm_bacterial_load$primary, Sm_bacterial_load$secondary_dose)
#THE GRAPH
#png(file="Sm_Pr.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Pr$day, Sm_Pr$status) ~ Sm_Pr$primary + Sm_Pr$secondary_dose), col=c("gray", "#D7EEF7", "#B5E8FC", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3,3,3,3, 3, 3, 3, 1,1,1,1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS","PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0","PBS-P.rett 3.0","PBS-P.rett 5.0","S.m-PBS", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#D7EEF7", "#B5E8FC", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Appendix Figure 1A - Controls
Sm_Pr_PBS<-subset(Sm_Pr,secondary_dose=="0")
#table(Sm_Pr_PBS$primary)
#THE GRAPH
#png(file="AFigure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=115)", "chronic infection only (n=114)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)
#dev.off()Figure 1 Panels
Figure 1A - 0.01
Sm_Pr_0.01<-subset(Sm_Pr,secondary_dose=="0.01")
#table(Sm_Pr_0.01$primary)
#THE GRAPH
#png(file="Figure1A.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=118)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1B - 0.1
Sm_Pr_0.1<-subset(Sm_Pr,secondary_dose=="0.1")
#table(Sm_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=206)", "chronic & secondary infection (n=196)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1c - 1.0
Sm_Pr_1.0<-subset(Sm_Pr,secondary_dose=="1")
#table(Sm_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure1C.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=283)", "chronic & secondary infection (n=275)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1D - 2.0
Sm_Pr_2.0<-subset(Sm_Pr,secondary_dose=="2")
#table(Sm_Pr_2.0$primary)
#THE GRAPH
#png(file="Figure1D.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_2.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=263)", "chronic & secondary infection (n=238)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1E - 3.0
Sm_Pr_3.0<-subset(Sm_Pr,secondary_dose=="3")
#table(Sm_Pr_3.0$primary)
#THE GRAPH
#png(file="Figure1E.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_3.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1F - 5.0
Sm_Pr_5.0<-subset(Sm_Pr,secondary_dose=="5")
#table(Sm_Pr_5.0$primary)
#THE GRAPH
#png(file="Figure1F.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Pr_5.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=165)", "chronic & secondary infection (n=159)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.7)
#dev.off()Figure 1
a <- ggdraw() +
draw_image("Figure1A.png")
b <- ggdraw() +
draw_image("Figure1B.png")
c <- ggdraw() +
draw_image("Figure1C.png")
d <- ggdraw() +
draw_image("Figure1D.png")
e <- ggdraw() +
draw_image("Figure1E.png")
f <- ggdraw() +
draw_image("Figure1F.png")
multiplot <- a + b + c + d + e + f
z <- multiplot +
plot_annotation(tag_levels = 'A')
ggsave("Fig.1.ChronicSm_Pr.pdf", z, device = "pdf")
z
Individual Blocks
Block 1
Sm_Pr_block1 <- subset(Sm_Pr, date == "8/2/2018")
table(Sm_Pr_block1$primary, Sm_Pr_block1$secondary_dose)##
## 0.01 0.1 1
## PBS 17 18 18
## S.m 15 15 15Sm_Pr_0.01_1<-subset(Sm_Pr_block1, secondary_dose=="0.01")
#table(Sm_Pr_0.01_1$primary)
Sm_Pr_0.1_1<-subset(Sm_Pr_block1, secondary_dose=="0.1")
#table(Sm_Pr_0.1_1$primary)
Sm_Pr_1.0_1<-subset(Sm_Pr_block1, secondary_dose=="1")
#table(Sm_Pr_1.0_1$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_3.0_1<-subset(Sm_Pr_block1, secondary_dose=="3")
#table(Sm_Pr_2.0_1$primary)
#Sm_Pr_5.0_1<-subset(Sm_Pr_block1, secondary_dose=="5")
#table(Sm_Pr_5.0_1$primary)
surv_Sm_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_1)
surv_Sm_Pr_0.01_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 8 5.3 1.37 3.19
## primary=S.m 15 3 5.7 1.28 3.19
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07surv_Sm_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_1)
surv_Sm_Pr_0.1_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18 18 15.5 0.395 2.33
## primary=S.m 15 12 14.5 0.423 2.33
##
## Chisq= 2.3 on 1 degrees of freedom, p= 0.1surv_Sm_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_1)
surv_Sm_Pr_1.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 18 18 15.8 0.301 5.3
## primary=S.m 15 13 15.2 0.314 5.3
##
## Chisq= 5.3 on 1 degrees of freedom, p= 0.02#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block1$day, Sm_Pr_block1$status) ~ Sm_Pr_block1$primary + Sm_Pr_block1$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Block 2
Sm_Pr_block2 <- subset(Sm_Pr, date == "8/6/2018")
table(Sm_Pr_block2$primary, Sm_Pr_block2$secondary_dose)##
## 0.01 0.1 1
## PBS 19 27 12
## S.m 24 23 15Sm_Pr_0.01_2<-subset(Sm_Pr_block2, secondary_dose=="0.01")
#table(Sm_Pr_0.01_2$primary)
Sm_Pr_0.1_2<-subset(Sm_Pr_block2, secondary_dose=="0.1")
#table(Sm_Pr_0.1_2$primary)
Sm_Pr_1.0_2<-subset(Sm_Pr_block2, secondary_dose=="1")
#table(Sm_Pr_1.0_2$primary)
#Sm_Pr_2.0_1<-subset(Sm_Pr_block1, secondary_dose=="2")
#table(Sm_Pr_2.0_1$primary)
surv_Sm_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_2)
surv_Sm_Pr_0.01_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 13 8.24 2.74 5.57
## primary=S.m 24 9 13.76 1.64 5.57
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02surv_Sm_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_2)
surv_Sm_Pr_0.1_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 20 13 3.72 10.5
## primary=S.m 23 8 15 3.24 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001surv_Sm_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_2)
surv_Sm_Pr_1.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 12 12 8.89 1.089 4.91
## primary=S.m 15 10 13.11 0.738 4.91
##
## Chisq= 4.9 on 1 degrees of freedom, p= 0.03#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block2$day, Sm_Pr_block2$status) ~ Sm_Pr_block2$primary + Sm_Pr_block2$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Block 3
Sm_Pr_block3 <- subset(Sm_Pr, date == "6/14/2019")
table(Sm_Pr_block3$primary, Sm_Pr_block3$secondary_dose)##
## 0.01 0.1 1 2
## PBS 15 14 17 20
## S.m 15 19 17 18Sm_Pr_0.01_3<-subset(Sm_Pr_block3, secondary_dose=="0.01")
#table(Sm_Pr_0.01_3$primary)
Sm_Pr_0.1_3<-subset(Sm_Pr_block3, secondary_dose=="0.1")
#table(Sm_Pr_0.1_3$primary)
Sm_Pr_1.0_3<-subset(Sm_Pr_block3, secondary_dose=="1")
#table(Sm_Pr_1.0_3$primary)
Sm_Pr_2.0_3<-subset(Sm_Pr_block3, secondary_dose=="2")
#table(Sm_Pr_2.0_3$primary)
surv_Sm_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_3)
surv_Sm_Pr_0.01_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 9 6.24 1.22 3.41
## primary=S.m 15 4 6.76 1.13 3.41
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06surv_Sm_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_3)
surv_Sm_Pr_0.1_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 14 7 2.74 6.63 12.1
## primary=S.m 19 0 4.26 4.26 12.1
##
## Chisq= 12.1 on 1 degrees of freedom, p= 5e-04surv_Sm_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_3)
surv_Sm_Pr_1.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 17 8.5 8.50 33
## primary=S.m 17 3 11.5 6.28 33
##
## Chisq= 33 on 1 degrees of freedom, p= 9e-09surv_Sm_Pr_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_3)
surv_Sm_Pr_2.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 12.1 5.15 26.8
## primary=S.m 18 5 12.9 4.83 26.8
##
## Chisq= 26.8 on 1 degrees of freedom, p= 2e-07#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block3$day, Sm_Pr_block3$status) ~ Sm_Pr_block3$primary + Sm_Pr_block3$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 4
Sm_Pr_block4 <- subset(Sm_Pr, date == "6/17/2019")
table(Sm_Pr_block4$primary, Sm_Pr_block4$secondary_dose)##
## 0.01 0.1 1 2
## PBS 15 19 17 17
## S.m 18 13 16 11Sm_Pr_0.01_4<-subset(Sm_Pr_block4, secondary_dose=="0.01")
#table(Sm_Pr_0.01_4$primary)
Sm_Pr_0.1_4<-subset(Sm_Pr_block4, secondary_dose=="0.1")
#table(Sm_Pr_0.1_4$primary)
Sm_Pr_1.0_4<-subset(Sm_Pr_block4, secondary_dose=="1")
#table(Sm_Pr_1.0_4$primary)
Sm_Pr_2.0_4<-subset(Sm_Pr_block4, secondary_dose=="2")
#table(Sm_Pr_2.0_4$primary)
surv_Sm_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_4)
surv_Sm_Pr_0.01_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 10 7.59 0.767 1.76
## primary=S.m 18 8 10.41 0.559 1.76
##
## Chisq= 1.8 on 1 degrees of freedom, p= 0.2surv_Sm_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_4)
surv_Sm_Pr_0.1_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 9 7.85 0.170 0.438
## primary=S.m 13 5 6.15 0.216 0.438
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5surv_Sm_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_4)
surv_Sm_Pr_1.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 16 9.76 4.00 15.7
## primary=S.m 16 6 12.24 3.18 15.7
##
## Chisq= 15.7 on 1 degrees of freedom, p= 7e-05surv_Sm_Pr_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_4)
surv_Sm_Pr_2.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 17 12.14 1.94 16.7
## primary=S.m 11 5 9.86 2.39 16.7
##
## Chisq= 16.7 on 1 degrees of freedom, p= 4e-05#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block4$day, Sm_Pr_block4$status) ~ Sm_Pr_block4$primary + Sm_Pr_block4$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 5
Sm_Pr_block5 <- subset(Sm_Pr, date == "6/19/2019")
table(Sm_Pr_block5$primary, Sm_Pr_block5$secondary_dose)##
## 0.01 0.1 1 2
## PBS 15 5 9 19
## S.m 11 17 14 15Sm_Pr_0.01_5<-subset(Sm_Pr_block5, secondary_dose=="0.01")
#table(Sm_Pr_0.01_5$primary)
Sm_Pr_0.1_5<-subset(Sm_Pr_block5, secondary_dose=="0.1")
#table(Sm_Pr_0.1_5$primary)
Sm_Pr_1.0_5<-subset(Sm_Pr_block5, secondary_dose=="1")
#table(Sm_Pr_1.0_5$primary)
Sm_Pr_2.0_5<-subset(Sm_Pr_block5, secondary_dose=="2")
#table(Sm_Pr_2.0_5$primary)
surv_Sm_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_5)
surv_Sm_Pr_0.01_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 10 7.83 0.604 1.61
## primary=S.m 11 5 7.17 0.659 1.61
##
## Chisq= 1.6 on 1 degrees of freedom, p= 0.2surv_Sm_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_5)
surv_Sm_Pr_0.1_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 5 5 1.21 11.90 17.1
## primary=S.m 17 5 8.79 1.64 17.1
##
## Chisq= 17.1 on 1 degrees of freedom, p= 4e-05surv_Sm_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_5)
surv_Sm_Pr_1.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 9 9 5.29 2.61 8.8
## primary=S.m 14 7 10.71 1.29 8.8
##
## Chisq= 8.8 on 1 degrees of freedom, p= 0.003surv_Sm_Pr_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_5)
surv_Sm_Pr_2.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 18 11.8 3.27 15
## primary=S.m 15 9 15.2 2.54 15
##
## Chisq= 15 on 1 degrees of freedom, p= 1e-04#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block5$day, Sm_Pr_block5$status) ~ Sm_Pr_block5$primary + Sm_Pr_block5$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 6
Sm_Pr_block6 <- subset(Sm_Pr, date == "6/25/2019")
table(Sm_Pr_block6$primary, Sm_Pr_block6$secondary_dose)##
## 0.01 0.1 1 2
## PBS 17 19 19 20
## S.m 20 15 17 16Sm_Pr_0.01_6<-subset(Sm_Pr_block6, secondary_dose=="0.01")
#table(Sm_Pr_0.01_6$primary)
Sm_Pr_0.1_6<-subset(Sm_Pr_block6, secondary_dose=="0.1")
#table(Sm_Pr_0.1_6$primary)
Sm_Pr_1.0_6<-subset(Sm_Pr_block6, secondary_dose=="1")
#table(Sm_Pr_1.0_6$primary)
Sm_Pr_2.0_6<-subset(Sm_Pr_block6, secondary_dose=="2")
#table(Sm_Pr_2.0_6$primary)
surv_Sm_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_6)
surv_Sm_Pr_0.01_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 12 7.11 3.36 7.48
## primary=S.m 20 5 9.89 2.42 7.48
##
## Chisq= 7.5 on 1 degrees of freedom, p= 0.006surv_Sm_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_6)
surv_Sm_Pr_0.1_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 15 8.52 4.92 13.2
## primary=S.m 15 2 8.48 4.95 13.2
##
## Chisq= 13.2 on 1 degrees of freedom, p= 3e-04surv_Sm_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_6)
surv_Sm_Pr_1.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 17 9.99 4.92 18
## primary=S.m 17 5 12.01 4.09 18
##
## Chisq= 18 on 1 degrees of freedom, p= 2e-05surv_Sm_Pr_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_6)
surv_Sm_Pr_2.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 11.9 5.57 27.8
## primary=S.m 16 6 14.1 4.68 27.8
##
## Chisq= 27.8 on 1 degrees of freedom, p= 1e-07#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block6$day, Sm_Pr_block6$status) ~ Sm_Pr_block6$primary + Sm_Pr_block6$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 7
Sm_Pr_block7 <- subset(Sm_Pr, date == "6/27/2019")
table(Sm_Pr_block7$primary, Sm_Pr_block7$secondary_dose)##
## 0.01 0.1 1 2
## PBS 20 17 20 19
## S.m 16 14 14 14Sm_Pr_0.01_7<-subset(Sm_Pr_block7, secondary_dose=="0.01")
#table(Sm_Pr_0.01_7$primary)
Sm_Pr_0.1_7<-subset(Sm_Pr_block7, secondary_dose=="0.1")
#table(Sm_Pr_0.1_7$primary)
Sm_Pr_1.0_7<-subset(Sm_Pr_block7, secondary_dose=="1")
#table(Sm_Pr_1.0_7$primary)
Sm_Pr_2.0_7<-subset(Sm_Pr_block7, secondary_dose=="2")
#table(Sm_Pr_2.0_7$primary)
surv_Sm_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.01_7)
surv_Sm_Pr_0.01_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 11 10.11 0.0779 0.206
## primary=S.m 16 9 9.89 0.0797 0.206
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7surv_Sm_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_7)
surv_Sm_Pr_0.1_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 17 14 9.38 2.27 6.94
## primary=S.m 14 6 10.62 2.01 6.94
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.008surv_Sm_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_7)
surv_Sm_Pr_1.0_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 13.5 3.09 22.5
## primary=S.m 14 7 13.5 3.11 22.5
##
## Chisq= 22.5 on 1 degrees of freedom, p= 2e-06surv_Sm_Pr_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_7)
surv_Sm_Pr_2.0_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 19 18 11.2 4.18 17.4
## primary=S.m 14 9 15.8 2.95 17.4
##
## Chisq= 17.4 on 1 degrees of freedom, p= 3e-05#surv_Sm_Pr_3.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_1)
#surv_Sm_Pr_3.0_1
#surv_Sm_Pr_5.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_1)
#surv_Sm_Pr_5.0_1
#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block7$day, Sm_Pr_block7$status) ~ Sm_Pr_block7$primary + Sm_Pr_block7$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "S.m-P.rett 0.01", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 8
Sm_Pr_block8 <- subset(Sm_Pr, date == "1/30/2020")
table(Sm_Pr_block8$primary, Sm_Pr_block8$secondary_dose)##
## 0 0.1 1 2 5
## PBS 18 29 27 26 28
## S.m 20 23 24 25 25Sm_Pr_0.1_8<-subset(Sm_Pr_block8, secondary_dose=="0.1")
#table(Sm_Pr_0.1_8$primary)
Sm_Pr_1.0_8<-subset(Sm_Pr_block8, secondary_dose=="1")
#table(Sm_Pr_1.0_8$primary)
Sm_Pr_2.0_8<-subset(Sm_Pr_block8, secondary_dose=="2")
#table(Sm_Pr_2.0_8$primary)
#Sm_Pr_3.0_8<-subset(Sm_Pr_block8, secondary_dose=="3")
#table(Sm_Pr_3.0_8$primary)
Sm_Pr_5.0_8<-subset(Sm_Pr_block8, secondary_dose=="5")
#table(Sm_Pr_5.0_8$primary)
surv_Sm_Pr_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_8)
surv_Sm_Pr_0.1_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 15.7 11.33 26.8
## primary=S.m 23 15 28.3 6.27 26.8
##
## Chisq= 26.8 on 1 degrees of freedom, p= 2e-07surv_Sm_Pr_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_8)
surv_Sm_Pr_1.0_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 16.9 5.97 33.4
## primary=S.m 24 19 29.1 3.48 33.4
##
## Chisq= 33.4 on 1 degrees of freedom, p= 8e-09surv_Sm_Pr_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_8)
surv_Sm_Pr_2.0_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26 26 15.3 7.49 36.4
## primary=S.m 25 20 30.7 3.73 36.4
##
## Chisq= 36.4 on 1 degrees of freedom, p= 2e-09#surv_Sm_Pr_3.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_8)
#surv_Sm_Pr_3.0_8
surv_Sm_Pr_5.0_8<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_8)
surv_Sm_Pr_5.0_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.3 10.49 48.2
## primary=S.m 25 16 28.7 5.61 48.2
##
## Chisq= 48.2 on 1 degrees of freedom, p= 4e-12#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block8$day, Sm_Pr_block8$status) ~ Sm_Pr_block8$primary + Sm_Pr_block8$secondary_dose), col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 9
Sm_Pr_block9 <- subset(Sm_Pr, date == "2/13/2020")
table(Sm_Pr_block9$primary, Sm_Pr_block9$secondary_dose)##
## 0 0.1 1 2 5
## PBS 18 28 31 30 30
## S.m 18 28 28 28 28Sm_Pr_0.1_9<-subset(Sm_Pr_block9, secondary_dose=="0.1")
#table(Sm_Pr_0.1_9$primary)
Sm_Pr_1.0_9<-subset(Sm_Pr_block9, secondary_dose=="1")
#table(Sm_Pr_1.0_9$primary)
Sm_Pr_2.0_9<-subset(Sm_Pr_block9, secondary_dose=="2")
#table(Sm_Pr_2.0_9$primary)
#Sm_Pr_3.0_9<-subset(Sm_Pr_block9, secondary_dose=="3")
#table(Sm_Pr_3.0_9$primary)
Sm_Pr_5.0_9<-subset(Sm_Pr_block9, secondary_dose=="5")
#table(Sm_Pr_5.0_9$primary)
surv_Sm_Pr_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_9)
surv_Sm_Pr_0.1_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 26 13.9 10.64 21.8
## primary=S.m 28 10 22.1 6.66 21.8
##
## Chisq= 21.8 on 1 degrees of freedom, p= 3e-06surv_Sm_Pr_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_9)
surv_Sm_Pr_1.0_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 31 31 18.9 7.72 41
## primary=S.m 28 15 27.1 5.39 41
##
## Chisq= 41 on 1 degrees of freedom, p= 2e-10surv_Sm_Pr_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_9)
surv_Sm_Pr_2.0_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 17.1 9.80 46.3
## primary=S.m 28 8 20.9 7.99 46.3
##
## Chisq= 46.3 on 1 degrees of freedom, p= 1e-11#surv_Sm_Pr_3.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_9)
#surv_Sm_Pr_3.0_9
surv_Sm_Pr_5.0_9<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_9)
surv_Sm_Pr_5.0_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 18.2 7.69 35.3
## primary=S.m 28 21 32.8 4.26 35.3
##
## Chisq= 35.3 on 1 degrees of freedom, p= 3e-09#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block9$day, Sm_Pr_block9$status) ~ Sm_Pr_block9$primary + Sm_Pr_block9$secondary_dose), col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 0.1", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 0.1", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 5.0"),col=c("gray", "#B5E8FC", "#09B6FC", "#055df5", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 10
Sm_Pr_block10 <- subset(Sm_Pr, date == "9/27/2020")
table(Sm_Pr_block10$primary, Sm_Pr_block10$secondary_dose)##
## 0 1 2 3 5
## PBS 10 27 26 20 20
## S.m 10 28 27 28 27#Sm_Pr_0.1_10<-subset(Sm_Pr_block10, secondary_dose=="0.1")
#table(Sm_Pr_0.1_10$primary)
Sm_Pr_1.0_10<-subset(Sm_Pr_block10, secondary_dose=="1")
#table(Sm_Pr_1.0_10$primary)
Sm_Pr_2.0_10<-subset(Sm_Pr_block10, secondary_dose=="2")
#table(Sm_Pr_2.0_10$primary)
Sm_Pr_3.0_10<-subset(Sm_Pr_block10, secondary_dose=="3")
#table(Sm_Pr_3.0_10$primary)
Sm_Pr_5.0_10<-subset(Sm_Pr_block10, secondary_dose=="5")
#table(Sm_Pr_5.0_10$primary)
#surv_Sm_Pr_0.1_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_10)
#surv_Sm_Pr_0.1_10
surv_Sm_Pr_1.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_10)
surv_Sm_Pr_1.0_10## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 12.7 16.16 53.6
## primary=S.m 28 7 21.3 9.62 53.6
##
## Chisq= 53.6 on 1 degrees of freedom, p= 2e-13surv_Sm_Pr_2.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_10)
surv_Sm_Pr_2.0_10## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 26 26 14.7 8.65 38.4
## primary=S.m 27 9 20.3 6.28 38.4
##
## Chisq= 38.4 on 1 degrees of freedom, p= 6e-10surv_Sm_Pr_3.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_10)
surv_Sm_Pr_3.0_10## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 10.4 8.82 30.9
## primary=S.m 28 10 19.6 4.69 30.9
##
## Chisq= 30.9 on 1 degrees of freedom, p= 3e-08surv_Sm_Pr_5.0_10<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_10)
surv_Sm_Pr_5.0_10## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_10)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 12.8 4.10 19.3
## primary=S.m 27 16 23.2 2.25 19.3
##
## Chisq= 19.3 on 1 degrees of freedom, p= 1e-05#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block10$day, Sm_Pr_block10$status) ~ Sm_Pr_block10$primary + Sm_Pr_block10$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 11
Sm_Pr_block11 <- subset(Sm_Pr, date == "10/6/2020")
table(Sm_Pr_block11$primary, Sm_Pr_block11$secondary_dose)##
## 0 1 2 3 5
## PBS 20 28 28 30 29
## S.m 16 28 27 26 27#Sm_Pr_0.1_11<-subset(Sm_Pr_block11, secondary_dose=="0.1")
#table(Sm_Pr_0.1_11$primary)
Sm_Pr_1.0_11<-subset(Sm_Pr_block11, secondary_dose=="1")
#table(Sm_Pr_1.0_11$primary)
Sm_Pr_2.0_11<-subset(Sm_Pr_block11, secondary_dose=="2")
#table(Sm_Pr_2.0_11$primary)
Sm_Pr_3.0_11<-subset(Sm_Pr_block11, secondary_dose=="3")
#table(Sm_Pr_3.0_11$primary)
Sm_Pr_5.0_11<-subset(Sm_Pr_block11, secondary_dose=="5")
#table(Sm_Pr_5.0_11$primary)
#surv_Sm_Pr_0.1_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_11)
#surv_Sm_Pr_0.1_11
surv_Sm_Pr_1.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_11)
surv_Sm_Pr_1.0_11## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 24 14.2 6.85 19
## primary=S.m 28 12 21.8 4.44 19
##
## Chisq= 19 on 1 degrees of freedom, p= 1e-05surv_Sm_Pr_2.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_11)
surv_Sm_Pr_2.0_11## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 14.4 12.93 51.5
## primary=S.m 27 9 22.6 8.21 51.5
##
## Chisq= 51.5 on 1 degrees of freedom, p= 7e-13surv_Sm_Pr_3.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_11)
surv_Sm_Pr_3.0_11## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 16.1 12.07 55
## primary=S.m 26 14 27.9 6.95 55
##
## Chisq= 55 on 1 degrees of freedom, p= 1e-13surv_Sm_Pr_5.0_11<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_11)
surv_Sm_Pr_5.0_11## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 16.1 10.4 47.6
## primary=S.m 27 11 23.9 7.0 47.6
##
## Chisq= 47.6 on 1 degrees of freedom, p= 5e-12#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block11$day, Sm_Pr_block11$status) ~ Sm_Pr_block11$primary + Sm_Pr_block11$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 12
Sm_Pr_block12 <- subset(Sm_Pr, date == "11/1/2020")
table(Sm_Pr_block12$primary, Sm_Pr_block12$secondary_dose)##
## 0 1 2 3 5
## PBS 20 30 30 30 30
## S.m 20 30 30 27 25#Sm_Pr_0.1_12<-subset(Sm_Pr_block12, secondary_dose=="0.1")
#table(Sm_Pr_0.1_12$primary)
Sm_Pr_1.0_12<-subset(Sm_Pr_block12, secondary_dose=="1")
#table(Sm_Pr_1.0_12$primary)
Sm_Pr_2.0_12<-subset(Sm_Pr_block12, secondary_dose=="2")
#table(Sm_Pr_2.0_12$primary)
Sm_Pr_3.0_12<-subset(Sm_Pr_block12, secondary_dose=="3")
#table(Sm_Pr_3.0_12$primary)
Sm_Pr_5.0_12<-subset(Sm_Pr_block12, secondary_dose=="5")
#table(Sm_Pr_5.0_12$primary)
#surv_Sm_Pr_0.1_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_12)
#surv_Sm_Pr_0.1_12
surv_Sm_Pr_1.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_12)
surv_Sm_Pr_1.0_12## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 28 13.6 15.17 33.2
## primary=S.m 30 13 27.4 7.55 33.2
##
## Chisq= 33.2 on 1 degrees of freedom, p= 8e-09surv_Sm_Pr_2.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_12)
surv_Sm_Pr_2.0_12## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 13.9 18.5 58.6
## primary=S.m 30 10 26.1 9.9 58.6
##
## Chisq= 58.6 on 1 degrees of freedom, p= 2e-14surv_Sm_Pr_3.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_12)
surv_Sm_Pr_3.0_12## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_3.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 15.6 13.39 50
## primary=S.m 27 21 35.4 5.88 50
##
## Chisq= 50 on 1 degrees of freedom, p= 2e-12surv_Sm_Pr_5.0_12<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_12)
surv_Sm_Pr_5.0_12## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_12)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 20.2 4.78 31.5
## primary=S.m 25 17 26.8 3.59 31.5
##
## Chisq= 31.5 on 1 degrees of freedom, p= 2e-08#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block12$day, Sm_Pr_block12$status) ~ Sm_Pr_block12$primary + Sm_Pr_block12$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 13
Sm_Pr_block13 <- subset(Sm_Pr, date == "6/16/2021")
table(Sm_Pr_block13$primary, Sm_Pr_block13$secondary_dose)##
## 0 0.1 1 2 5
## PBS 29 30 28 28 28
## S.m 30 29 29 27 27Sm_Pr_0.1_13<-subset(Sm_Pr_block13, secondary_dose=="0.1")
#table(Sm_Pr_0.1_13$primary)
Sm_Pr_1.0_13<-subset(Sm_Pr_block13, secondary_dose=="1")
#table(Sm_Pr_1.0_13$primary)
Sm_Pr_2.0_13<-subset(Sm_Pr_block13, secondary_dose=="2")
#table(Sm_Pr_2.0_13$primary)
#Sm_Pr_3.0_13<-subset(Sm_Pr_block13, secondary_dose=="3")
#table(Sm_Pr_2.0_13$primary)
Sm_Pr_5.0_13<-subset(Sm_Pr_block13, secondary_dose=="5")
#table(Sm_Pr_2.0_13$primary)
surv_Sm_Pr_0.1_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_0.1_13)
surv_Sm_Pr_0.1_13## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_0.1_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 28 14.5 12.56 39.3
## primary=S.m 29 6 19.5 9.34 39.3
##
## Chisq= 39.3 on 1 degrees of freedom, p= 4e-10surv_Sm_Pr_1.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_1.0_13)
surv_Sm_Pr_1.0_13## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_1.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 14.3 13.2 51.4
## primary=S.m 29 4 17.7 10.6 51.4
##
## Chisq= 51.4 on 1 degrees of freedom, p= 7e-13surv_Sm_Pr_2.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_2.0_13)
surv_Sm_Pr_2.0_13## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_2.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.4 10.42 45.1
## primary=S.m 27 5 17.6 9.06 45.1
##
## Chisq= 45.1 on 1 degrees of freedom, p= 2e-11#surv_Sm_Pr_3.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_3.0_13)
#surv_Sm_Pr_3.0_13
surv_Sm_Pr_5.0_13<-survdiff(Surv(day, status) ~ primary, data=Sm_Pr_5.0_13)
surv_Sm_Pr_5.0_13## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Pr_5.0_13)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 15.4 10.42 45.1
## primary=S.m 27 4 16.6 9.61 45.1
##
## Chisq= 45.1 on 1 degrees of freedom, p= 2e-11#THE GRAPH
par(mar=c(4,5,1,1))
plot(survfit(Surv(Sm_Pr_block13$day, Sm_Pr_block13$status) ~ Sm_Pr_block13$primary + Sm_Pr_block13$secondary_dose), col=c("gray", "#09B6FC", "#055df5", "#003694", "black"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("PBS-PBS", "PBS-P.rett 1.0", "PBS-P.rett 2.0", "PBS-P.rett 3.0", "PBS-P.rett 5.0", "S.m-PBS", "S.m-P.rett 1.0", "S.m-P.rett 2.0", "S.m-P.rett 3.0", "S.m-P.rett 5.0"),col=c("gray", "#09B6FC", "#055df5", "#003694", "black"),lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=5, cex=1)
Statistical tests
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_0.01_1<-pchisq(surv_Sm_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_0.1_1<-pchisq(surv_Sm_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_1.0_1<-pchisq(surv_Sm_Pr_1.0_1$chisq, df=1, lower.tail=F)
pvalue_0.01_2<-pchisq(surv_Sm_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_0.1_2<-pchisq(surv_Sm_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_1.0_2<-pchisq(surv_Sm_Pr_1.0_2$chisq, df=1, lower.tail=F)
pvalue_0.01_3<-pchisq(surv_Sm_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_0.1_3<-pchisq(surv_Sm_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_1.0_3<-pchisq(surv_Sm_Pr_1.0_3$chisq, df=1, lower.tail=F)
pvalue_2.0_3<-pchisq(surv_Sm_Pr_2.0_3$chisq, df=1, lower.tail=F)
pvalue_0.01_4<-pchisq(surv_Sm_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_0.1_4<-pchisq(surv_Sm_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_1.0_4<-pchisq(surv_Sm_Pr_1.0_4$chisq, df=1, lower.tail=F)
pvalue_2.0_4<-pchisq(surv_Sm_Pr_2.0_4$chisq, df=1, lower.tail=F)
pvalue_0.01_5<-pchisq(surv_Sm_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_0.1_5<-pchisq(surv_Sm_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_1.0_5<-pchisq(surv_Sm_Pr_1.0_5$chisq, df=1, lower.tail=F)
pvalue_2.0_5<-pchisq(surv_Sm_Pr_2.0_5$chisq, df=1, lower.tail=F)
pvalue_0.01_6<-pchisq(surv_Sm_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_0.1_6<-pchisq(surv_Sm_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_1.0_6<-pchisq(surv_Sm_Pr_1.0_6$chisq, df=1, lower.tail=F)
pvalue_2.0_6<-pchisq(surv_Sm_Pr_2.0_6$chisq, df=1, lower.tail=F)
pvalue_0.01_7<-pchisq(surv_Sm_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_0.1_7<-pchisq(surv_Sm_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_1.0_7<-pchisq(surv_Sm_Pr_1.0_7$chisq, df=1, lower.tail=F)
pvalue_2.0_7<-pchisq(surv_Sm_Pr_2.0_7$chisq, df=1, lower.tail=F)
pvalue_0.1_8<-pchisq(surv_Sm_Pr_0.1_8$chisq, df=1, lower.tail=F)
pvalue_1.0_8<-pchisq(surv_Sm_Pr_1.0_8$chisq, df=1, lower.tail=F)
pvalue_2.0_8<-pchisq(surv_Sm_Pr_2.0_8$chisq, df=1, lower.tail=F)
pvalue_5.0_8<-pchisq(surv_Sm_Pr_5.0_8$chisq, df=1, lower.tail=F)
pvalue_0.1_9<-pchisq(surv_Sm_Pr_0.1_9$chisq, df=1, lower.tail=F)
pvalue_1.0_9<-pchisq(surv_Sm_Pr_1.0_9$chisq, df=1, lower.tail=F)
pvalue_2.0_9<-pchisq(surv_Sm_Pr_2.0_9$chisq, df=1, lower.tail=F)
pvalue_5.0_9<-pchisq(surv_Sm_Pr_5.0_9$chisq, df=1, lower.tail=F)
pvalue_1.0_10<-pchisq(surv_Sm_Pr_1.0_10$chisq, df=1, lower.tail=F)
pvalue_2.0_10<-pchisq(surv_Sm_Pr_2.0_10$chisq, df=1, lower.tail=F)
pvalue_3.0_10<-pchisq(surv_Sm_Pr_3.0_10$chisq, df=1, lower.tail=F)
pvalue_5.0_10<-pchisq(surv_Sm_Pr_5.0_10$chisq, df=1, lower.tail=F)
pvalue_1.0_11<-pchisq(surv_Sm_Pr_1.0_11$chisq, df=1, lower.tail=F)
pvalue_2.0_11<-pchisq(surv_Sm_Pr_2.0_11$chisq, df=1, lower.tail=F)
pvalue_3.0_11<-pchisq(surv_Sm_Pr_3.0_11$chisq, df=1, lower.tail=F)
pvalue_5.0_11<-pchisq(surv_Sm_Pr_5.0_11$chisq, df=1, lower.tail=F)
pvalue_1.0_12<-pchisq(surv_Sm_Pr_1.0_12$chisq, df=1, lower.tail=F)
pvalue_2.0_12<-pchisq(surv_Sm_Pr_2.0_12$chisq, df=1, lower.tail=F)
pvalue_3.0_12<-pchisq(surv_Sm_Pr_3.0_12$chisq, df=1, lower.tail=F)
pvalue_5.0_12<-pchisq(surv_Sm_Pr_5.0_12$chisq, df=1, lower.tail=F)
pvalue_0.1_13<-pchisq(surv_Sm_Pr_0.1_13$chisq, df=1, lower.tail=F)
pvalue_1.0_13<-pchisq(surv_Sm_Pr_1.0_13$chisq, df=1, lower.tail=F)
pvalue_2.0_13<-pchisq(surv_Sm_Pr_2.0_13$chisq, df=1, lower.tail=F)
pvalue_5.0_13<-pchisq(surv_Sm_Pr_5.0_13$chisq, df=1, lower.tail=F)# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_0.01 <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7)
test.stat_0.01 <- -2*sum(log(pvalues_0.01))
pfinal_0.01<-pchisq(test.stat_0.01, df=2*length(pvalues_0.01), lower.tail=F)
pvalues_0.1 <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9, pvalue_0.1_13)
test.stat_0.1 <- -2*sum(log(pvalues_0.1))
pfinal_0.1<-pchisq(test.stat_0.1, df=2*length(pvalues_0.1), lower.tail=F)
pvalues_1.0 <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
test.stat_1.0 <- -2*sum(log(pvalues_1.0))
pfinal_1.0<-pchisq(test.stat_1.0, df=2*length(pvalues_1.0), lower.tail=F)
pvalues_2.0 <- c(pvalue_2.0_3,pvalue_2.0_4,pvalue_2.0_5,pvalue_2.0_6,pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
test.stat_2.0 <- -2*sum(log(pvalues_2.0))
pfinal_2.0<-pchisq(test.stat_2.0, df=2*length(pvalues_2.0), lower.tail=F)
pvalues_3.0 <- c(pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12)
test.stat_3.0 <- -2*sum(log(pvalues_3.0))
pfinal_3.0<-pchisq(test.stat_3.0, df=2*length(pvalues_3.0), lower.tail=F)
pvalues_5.0 <- c(pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
test.stat_5.0 <- -2*sum(log(pvalues_5.0))
pfinal_5.0<-pchisq(test.stat_5.0, df=2*length(pvalues_5.0), lower.tail=F)# Generate Table of final results
pvalues_0.01_table <- c(pvalue_0.01_1,pvalue_0.01_2,pvalue_0.01_3,pvalue_0.01_4,pvalue_0.01_5,pvalue_0.01_6,pvalue_0.01_7, NA, NA, NA, NA, NA, NA)
p0.01<-scales::pvalue(pvalues_0.01_table)
p0.01 <- p0.01 %>% replace_na("--")
pvalues_0.1_table <- c(pvalue_0.1_1,pvalue_0.1_2,pvalue_0.1_3,pvalue_0.1_4,pvalue_0.1_5,pvalue_0.1_6,pvalue_0.1_7, pvalue_0.1_8, pvalue_0.1_9,NA, NA, NA, pvalue_0.1_13)
p0.1<-scales::pvalue(pvalues_0.1_table)
p0.1<- p0.1 %>% replace_na("--")
pvalues_1.0_table <- c(pvalue_1.0_1,pvalue_1.0_2,pvalue_1.0_3,pvalue_1.0_4,pvalue_1.0_5,pvalue_1.0_6,pvalue_1.0_7, pvalue_1.0_8, pvalue_1.0_9, pvalue_1.0_10, pvalue_1.0_11, pvalue_1.0_12, pvalue_1.0_13)
p1.0<-scales::pvalue(pvalues_1.0_table)
pvalues_2.0_table <- c(NA, NA, pvalue_2.0_3, pvalue_2.0_4, pvalue_2.0_5, pvalue_2.0_6, pvalue_2.0_7, pvalue_2.0_8, pvalue_2.0_9, pvalue_2.0_10, pvalue_2.0_11, pvalue_2.0_12, pvalue_2.0_13)
p2.0<-scales::pvalue(pvalues_2.0_table)
p2.0<-p2.0 %>% replace_na("--")
pvalues_3.0_table <- c(NA, NA, NA, NA, NA, NA, NA, NA, NA, pvalue_3.0_10, pvalue_3.0_11, pvalue_3.0_12, NA)
p3.0<-scales::pvalue(pvalues_3.0_table)
p3.0<-p3.0 %>% replace_na("--")
pvalues_5.0_table <- c(NA, NA, NA, NA, NA, NA, NA, pvalue_5.0_8, pvalue_5.0_9, pvalue_5.0_10, pvalue_5.0_11, pvalue_5.0_12, pvalue_5.0_13)
p5.0<-scales::pvalue(pvalues_5.0_table)
p5.0<-p5.0 %>% replace_na("--")
pvalue_table <- data.frame("0.01" = p0.01, "0.1" = p0.1, "1.0" = p1.0, "2.0" = p2.0, "3.0" = p3.0, "5.0" = p5.0 )
colnames(pvalue_table) <- c(0.01, 0.1, 1.0, 2.0, 3.0, 5.0)
pfisher<- c(pfinal_0.01, pfinal_0.1, pfinal_1.0, pfinal_2.0, pfinal_3.0, pfinal_5.0)
pfisher<-scales::pvalue(pfisher)
pvalue_table[nrow(pvalue_table) + 1,]<-pfisher
row.names(pvalue_table) <- c("Block 1 - August 2nd, 2018", "Block 2 - August 6th, 2018", "Block 3 - June 14, 2019", "Block 4 - June 17, 2019", "Block 5 - June 19, 2019", "Block 6 - June 25, 2019", "Block 7 - June 27, 2019", "Block 8 - January 30, 2020", "Block 9 - February 13th, 2020", "Block 10 - September 27th, 2020", "Block 11 - October 6th, 2020", "Block 12 - November 1st, 2020", "Block 13 - June 16th, 2021", " Final P-value ")
pvalue_table %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 6)) %>%
pack_rows("Individual Blocks", 1, 13) %>%
pack_rows("Fisher's Combined", 14, 14)|
Infectious Dose (Abs600nm)
|
||||||
|---|---|---|---|---|---|---|
| 0.01 | 0.1 | 1 | 2 | 3 | 5 | |
| Individual Blocks | ||||||
| Block 1 - August 2nd, 2018 | 0.074 | 0.127 | 0.021 | – | – | – |
| Block 2 - August 6th, 2018 | 0.018 | 0.001 | 0.027 | – | – | – |
| Block 3 - June 14, 2019 | 0.065 | <0.001 | <0.001 | <0.001 | – | – |
| Block 4 - June 17, 2019 | 0.185 | 0.508 | <0.001 | <0.001 | – | – |
| Block 5 - June 19, 2019 | 0.205 | <0.001 | 0.003 | <0.001 | – | – |
| Block 6 - June 25, 2019 | 0.006 | <0.001 | <0.001 | <0.001 | – | – |
| Block 7 - June 27, 2019 | 0.650 | 0.008 | <0.001 | <0.001 | – | – |
| Block 8 - January 30, 2020 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Block 9 - February 13th, 2020 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Block 10 - September 27th, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 11 - October 6th, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 12 - November 1st, 2020 | – | – | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 13 - June 16th, 2021 | – | <0.001 | <0.001 | <0.001 | – | <0.001 |
| Fisher’s Combined | ||||||
| Final P-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 |
- Conclusion: There was significant protection at all infectious doses for the secondary infection (p<0.001)
Figure 2 - E. faecalis Chronic Infection
- Mortality over 7 days following secondary infection with P. rettgeri in flies chronically infected with E. faecalis versus control flies not given a chronic infection.
Graph
#KAPLAN MEIER
#table(Ef_Pr$primary, Ef_Pr$secondary_dose) #tells you the number of flies in each condition
#png(file="E.fae-P.rett.png", width=1200, heigh=900, res=130) #determines where the graph will be saved
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
#The graph
plot(survfit(Surv(Ef_Pr$day, Ef_Pr$status) ~ Ef_Pr$primary + Ef_Pr$secondary_dose),
col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #line color, cycles through
lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type
lwd=5, #line width
yaxt='n', #hide the y axis so I can customize it
xaxt='n', #hide the x axis so I can customize itS
ylab="proportion alive", #label for the y-axis
xlab="days post-secondary infection", #label for the x-axis
cex.lab=2.5, #font size of axis labels
xlim=c(0,7), #determines range of the x-axis
xaxs='i') #makes the data flush with the axis
#Y axis
axis(2, #left side
las=2, #perpendicular to the axis
cex.axis=1.5, #font size
lwd=5) #line width
#X axis
axis(1, #bottom of the graph
cex.axis=1.5, #font size
lwd=5) #line width
box(lwd=5)
#The legend
legend("topright", #location
bty="n", #no box around the legend
c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"), #labels in the legend
col=c("#B5E8FC", "#09B6FC", "#055df5", "#0600a8"), #colors for the legend lines, cycles through
lty=c(3, 3, 3, 3, 1, 1, 1, 1), #line type for legend lines
lwd=5, #line width for legend lines
cex=1) #font size
Figure 2 Panels
Figure 2A - 0.001
Ef_Pr_0.001<-subset(Ef_Pr,secondary_dose=="0.001")
#table(Ef_Pr_0.001$primary)
#THE GRAPH
#png(file="Figure2A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=101)", "chronic & secondary infection (n=110)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()Figure 2B - 0.01
Ef_Pr_0.01<-subset(Ef_Pr,secondary_dose=="0.01")
#table(Ef_Pr_0.01$primary)
#THE GRAPH
#png(file="Figure2B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=148)", "chronic & secondary infection (n=119)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()Figure 2C - 0.1
Ef_Pr_0.1<-subset(Ef_Pr,secondary_dose=="0.1")
#table(Ef_Pr_0.1$primary)
#THE GRAPH
#png(file="Figure2C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=160)", "chronic & secondary infection (n=115)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()Figure 2D - 1.0
Ef_Pr_1.0<-subset(Ef_Pr,secondary_dose=="1")
#table(Ef_Pr_1.0$primary)
#THE GRAPH
#png(file="Figure2D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Ef_Pr_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=157)", "chronic & secondary infection (n=125)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.8)
#dev.off()Figure 2
g <- ggdraw() +
draw_image("Figure2A.png")
h <- ggdraw() +
draw_image("Figure2B.png")
i <- ggdraw() +
draw_image("Figure2C.png")
j <- ggdraw() +
draw_image("Figure2D.png")
multiplot <- g + h + i + j
y <- multiplot +
plot_annotation(tag_levels = 'A')
ggsave("Fig.2.ChronicEf_Pr.pdf", y, device = "pdf")
y
Individual Blocks
Block 1
Ef_Pr_block1 <- subset(Ef_Pr, Date == "6/11/2019")
table(Ef_Pr_block1 $primary, Ef_Pr_block1 $secondary_dose)##
## 0.01 0.1 1
## 01_PBS 20 19 20
## 02_Ef 16 16 17#Ef_Pr_0.001_1<-subset(Ef_Pr_block1, secondary_dose=="2")
#table(Ef_Pr_0.001_1$primary)
Ef_Pr_0.01_1<-subset(Ef_Pr_block1, secondary_dose=="0.01")
#table(Ef_Pr_0.01_1$primary)
Ef_Pr_0.1_1<-subset(Ef_Pr_block1, secondary_dose=="0.1")
#table(Ef_Pr_0.1_1$primary)
Ef_Pr_1.0_1<-subset(Ef_Pr_block1, secondary_dose=="1")
#table(Ef_Pr_1.0_1$primary)
surv_Ef_Pr_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_1)
surv_Ef_Pr_0.01_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 19 12.8 3.03 8.21
## primary=02_Ef 16 16 22.2 1.74 8.21
##
## Chisq= 8.2 on 1 degrees of freedom, p= 0.004surv_Ef_Pr_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_1)
surv_Ef_Pr_0.1_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19 19 14.2 1.66 5.9
## primary=02_Ef 16 15 19.8 1.18 5.9
##
## Chisq= 5.9 on 1 degrees of freedom, p= 0.02surv_Ef_Pr_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_1)
surv_Ef_Pr_1.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 20 17.8 0.262 5.13
## primary=02_Ef 17 17 19.2 0.244 5.13
##
## Chisq= 5.1 on 1 degrees of freedom, p= 0.02par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block1$day, Ef_Pr_block1$status) ~ Ef_Pr_block1$primary + Ef_Pr_block1$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Block 2
Ef_Pr_block2 <- subset(Ef_Pr, Date == "6/22/2018")
table(Ef_Pr_block2 $primary, Ef_Pr_block2 $secondary_dose)##
## 0.001 0.01 0.1 1
## 01_PBS 28 18 29 30
## 02_Ef 33 18 19 29Ef_Pr_0.001_2<-subset(Ef_Pr_block2, secondary_dose=="0.001")
table(Ef_Pr_0.001_2$primary)##
## 01_PBS 02_Ef
## 28 33Ef_Pr_0.01_2<-subset(Ef_Pr_block2, secondary_dose=="0.01")
#table(Ef_Pr_0.01_2$primary)
Ef_Pr_0.1_2<-subset(Ef_Pr_block2, secondary_dose=="0.1")
#table(Ef_Pr_0.1_2$primary)
Ef_Pr_1.0_2<-subset(Ef_Pr_block2, secondary_dose=="1")
#table(Ef_Pr_1.0_2$primary)
surv_Ef_Pr_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_2)
surv_Ef_Pr_0.001_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 28 13 10.9 0.406 1.01
## primary=02_Ef 33 11 13.1 0.337 1.01
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3surv_Ef_Pr_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_2)
surv_Ef_Pr_0.01_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 13 8.44 2.46 6.03
## primary=02_Ef 18 6 10.56 1.97 6.03
##
## Chisq= 6 on 1 degrees of freedom, p= 0.01surv_Ef_Pr_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_2)
surv_Ef_Pr_0.1_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 23 17.1 2.01 7.88
## primary=02_Ef 19 8 13.9 2.48 7.88
##
## Chisq= 7.9 on 1 degrees of freedom, p= 0.005surv_Ef_Pr_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_2)
surv_Ef_Pr_1.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 30 30 27.5 0.235 5.56
## primary=02_Ef 29 26 28.5 0.226 5.56
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block2$day, Ef_Pr_block2$status) ~ Ef_Pr_block2$primary + Ef_Pr_block2$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 3
Ef_Pr_block3 <- subset(Ef_Pr, Date == "6/25/2018")
table(Ef_Pr_block3 $primary, Ef_Pr_block3 $secondary_dose)##
## 0.001 0.01 0.1 1
## 01_PBS 25 27 29 37
## 02_Ef 23 15 25 19Ef_Pr_0.001_3<-subset(Ef_Pr_block3, secondary_dose=="0.001")
table(Ef_Pr_0.001_3$primary)##
## 01_PBS 02_Ef
## 25 23Ef_Pr_0.01_3<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_3$primary)
Ef_Pr_0.1_3<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_3$primary)
Ef_Pr_1.0_3<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_3$primary)
surv_Ef_Pr_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_3)
surv_Ef_Pr_0.001_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 25 18 14.4 0.912 3.13
## primary=02_Ef 23 11 14.6 0.897 3.13
##
## Chisq= 3.1 on 1 degrees of freedom, p= 0.08surv_Ef_Pr_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_3)
surv_Ef_Pr_0.01_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27 21 15.8 1.68 7.31
## primary=02_Ef 15 5 10.2 2.62 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007surv_Ef_Pr_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_3)
surv_Ef_Pr_0.1_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 21 20.6 0.00806 0.0384
## primary=02_Ef 25 18 18.4 0.00902 0.0384
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8surv_Ef_Pr_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_3)
surv_Ef_Pr_1.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37 36 35.7 0.00290 0.235
## primary=02_Ef 19 18 18.3 0.00564 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block3$day, Ef_Pr_block3$status) ~ Ef_Pr_block3$primary + Ef_Pr_block3$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 4
Not sure why the data stops after 3 days, Frank did this one
Ef_Pr_block4 <- subset(Ef_Pr, Date == "7/19/2018")
table(Ef_Pr_block4 $primary, Ef_Pr_block4 $secondary_dose)##
## 0.001 0.01 0.1 1
## 01_PBS 9 18 18 10
## 02_Ef 23 27 16 17Ef_Pr_0.001_4<-subset(Ef_Pr_block4, secondary_dose=="0.001")
#table(Ef_Pr_0.001_4$primary)
Ef_Pr_0.01_4<-subset(Ef_Pr_block4, secondary_dose=="0.01")
#table(Ef_Pr_0.01_4$primary)
Ef_Pr_0.1_4<-subset(Ef_Pr_block4, secondary_dose=="0.1")
#table(Ef_Pr_0.1_4$primary)
Ef_Pr_1.0_4<-subset(Ef_Pr_block4, secondary_dose=="1")
#table(Ef_Pr_1.0_4$primary)
surv_Ef_Pr_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_4)
surv_Ef_Pr_0.001_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 9 3 2.56 0.0743 0.123
## primary=02_Ef 23 6 6.44 0.0296 0.123
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7surv_Ef_Pr_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_4)
surv_Ef_Pr_0.01_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 12 12 7.93e-05 0.00026
## primary=02_Ef 27 16 16 5.98e-05 0.00026
##
## Chisq= 0 on 1 degrees of freedom, p= 1surv_Ef_Pr_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_4)
surv_Ef_Pr_0.1_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 9 8.31 0.0574 0.147
## primary=02_Ef 16 9 9.69 0.0492 0.147
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7surv_Ef_Pr_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_4)
surv_Ef_Pr_1.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 10 9 8.79 0.00516 0.0292
## primary=02_Ef 17 17 17.21 0.00263 0.0292
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block4$day, Ef_Pr_block4$status) ~ Ef_Pr_block4$primary + Ef_Pr_block4$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 5
Ef_Pr_block5 <- subset(Ef_Pr, Date == "6/3/2019")
table(Ef_Pr_block5 $primary, Ef_Pr_block5 $secondary_dose)##
## 0.001 0.01 0.1 1
## 01_PBS 18 21 19 17
## 02_Ef 12 13 6 12Ef_Pr_0.001_5<-subset(Ef_Pr_block5, secondary_dose=="0.001")
table(Ef_Pr_0.001_5$primary)##
## 01_PBS 02_Ef
## 18 12Ef_Pr_0.01_5<-subset(Ef_Pr_block5, secondary_dose=="0.01")
#table(Ef_Pr_0.01_5$primary)
Ef_Pr_0.1_5<-subset(Ef_Pr_block5, secondary_dose=="0.1")
#table(Ef_Pr_0.1_5$primary)
Ef_Pr_1.0_5<-subset(Ef_Pr_block5, secondary_dose=="1")
#table(Ef_Pr_1.0_5$primary)
surv_Ef_Pr_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_5)
surv_Ef_Pr_0.001_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 18 13 10.65 0.518 1.39
## primary=02_Ef 12 6 8.35 0.661 1.39
##
## Chisq= 1.4 on 1 degrees of freedom, p= 0.2surv_Ef_Pr_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_5)
surv_Ef_Pr_0.01_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 16 12.9 0.759 2.31
## primary=02_Ef 13 7 10.1 0.965 2.31
##
## Chisq= 2.3 on 1 degrees of freedom, p= 0.1surv_Ef_Pr_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_5)
surv_Ef_Pr_0.1_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 19 15 14.05 0.0636 0.502
## primary=02_Ef 6 4 4.95 0.1806 0.502
##
## Chisq= 0.5 on 1 degrees of freedom, p= 0.5surv_Ef_Pr_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_5)
surv_Ef_Pr_1.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 17 17 14.7 0.375 6.35
## primary=02_Ef 12 12 14.3 0.383 6.35
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block5$day, Ef_Pr_block5$status) ~ Ef_Pr_block5$primary + Ef_Pr_block5$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 6
Ef_Pr_block6 <- subset(Ef_Pr, Date == "6/5/2019")
table(Ef_Pr_block6 $primary, Ef_Pr_block6 $secondary_dose)##
## 0.001 0.01 0.1 1
## 01_PBS 21 21 22 20
## 02_Ef 19 16 16 17Ef_Pr_0.001_6<-subset(Ef_Pr_block6, secondary_dose=="0.001")
table(Ef_Pr_0.001_6$primary)##
## 01_PBS 02_Ef
## 21 19Ef_Pr_0.01_6<-subset(Ef_Pr_block6, secondary_dose=="0.01")
#table(Ef_Pr_0.01_6$primary)
Ef_Pr_0.1_6<-subset(Ef_Pr_block6, secondary_dose=="0.1")
#table(Ef_Pr_0.1_6$primary)
Ef_Pr_1.0_6<-subset(Ef_Pr_block6, secondary_dose=="1")
#table(Ef_Pr_1.0_6$primary)
surv_Ef_Pr_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_6)
surv_Ef_Pr_0.001_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 13 8.18 2.84 6.64
## primary=02_Ef 19 4 8.82 2.63 6.64
##
## Chisq= 6.6 on 1 degrees of freedom, p= 0.01surv_Ef_Pr_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_6)
surv_Ef_Pr_0.01_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 21 10 7.63 0.738 1.76
## primary=02_Ef 16 4 6.37 0.883 1.76
##
## Chisq= 1.8 on 1 degrees of freedom, p= 0.2surv_Ef_Pr_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_6)
surv_Ef_Pr_0.1_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 22 10 10.96 0.0834 0.246
## primary=02_Ef 16 9 8.04 0.1137 0.246
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6surv_Ef_Pr_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_6)
surv_Ef_Pr_1.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 20 17 16.7 0.00591 0.026
## primary=02_Ef 17 16 16.3 0.00604 0.026
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block6$day, Ef_Pr_block6$status) ~ Ef_Pr_block6$primary + Ef_Pr_block6$secondary_dose), col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 3, 1, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.001", "PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.001", "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#B5E8FC", "#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 3, 1, 1, 1, 1),lwd=5, cex=1)
Block 7
Ef_Pr_block7 <- subset(Ef_Pr, Date == "7/1/2019")
table(Ef_Pr_block7 $primary, Ef_Pr_block7 $secondary_dose)##
## 0.01 0.1 1
## 01_PBS 23 24 23
## 02_Ef 14 17 14#Ef_Pr_0.001_7<-subset(Ef_Pr_block3, secondary_dose=="0.001")
#table(Ef_Pr_0.001_7$primary)
Ef_Pr_0.01_7<-subset(Ef_Pr_block3, secondary_dose=="0.01")
#table(Ef_Pr_0.01_7$primary)
Ef_Pr_0.1_7<-subset(Ef_Pr_block3, secondary_dose=="0.1")
#table(Ef_Pr_0.1_7$primary)
Ef_Pr_1.0_7<-subset(Ef_Pr_block3, secondary_dose=="1")
#table(Ef_Pr_1.0_7$primary)
#surv_Ef_Pr_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.001_7)
#surv_Ef_Pr_0.001_7
surv_Ef_Pr_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.01_7)
surv_Ef_Pr_0.01_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 27 21 15.8 1.68 7.31
## primary=02_Ef 15 5 10.2 2.62 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007surv_Ef_Pr_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_0.1_7)
surv_Ef_Pr_0.1_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 29 21 20.6 0.00806 0.0384
## primary=02_Ef 25 18 18.4 0.00902 0.0384
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8surv_Ef_Pr_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Ef_Pr_1.0_7)
surv_Ef_Pr_1.0_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ef_Pr_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=01_PBS 37 36 35.7 0.00290 0.235
## primary=02_Ef 19 18 18.3 0.00564 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6par(mar=c(4,5,1,1))
plot(survfit(Surv(Ef_Pr_block7$day, Ef_Pr_block7$status) ~ Ef_Pr_block7$primary + Ef_Pr_block7$secondary_dose), col=c("#09B6FC", "#055df5", "#003694"), lty=c(3, 3, 3, 1, 1, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("topright",inset=0, bty="n" , c("PBS-P.rett 0.01", "PBS-P.rett 0.1", 'PBS-P.rett 1.0', "E.fae-P.rett 0.01", "E.fae-P.rett 0.1", "E.fae-P.rett 1.0"),col=c("#09B6FC", "#055df5", "#003694"),lty=c(3, 3, 3, 1, 1, 1),lwd=5, cex=1)
Statistical Tests
- In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ef_Pr_0.01_1<-pchisq(surv_Ef_Pr_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_1<-pchisq(surv_Ef_Pr_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_1<-pchisq(surv_Ef_Pr_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_2<-pchisq(surv_Ef_Pr_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_2<-pchisq(surv_Ef_Pr_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_2<-pchisq(surv_Ef_Pr_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_2<-pchisq(surv_Ef_Pr_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_3<-pchisq(surv_Ef_Pr_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_3<-pchisq(surv_Ef_Pr_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_3<-pchisq(surv_Ef_Pr_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_3<-pchisq(surv_Ef_Pr_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_4<-pchisq(surv_Ef_Pr_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_4<-pchisq(surv_Ef_Pr_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_4<-pchisq(surv_Ef_Pr_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_4<-pchisq(surv_Ef_Pr_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_5<-pchisq(surv_Ef_Pr_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_5<-pchisq(surv_Ef_Pr_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_5<-pchisq(surv_Ef_Pr_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_5<-pchisq(surv_Ef_Pr_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.001_6<-pchisq(surv_Ef_Pr_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_6<-pchisq(surv_Ef_Pr_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_6<-pchisq(surv_Ef_Pr_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_6<-pchisq(surv_Ef_Pr_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.01_7<-pchisq(surv_Ef_Pr_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_0.1_7<-pchisq(surv_Ef_Pr_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Ef_Pr_1.0_7<-pchisq(surv_Ef_Pr_1.0_7$chisq, df=1, lower.tail=F)# Use all p-values from a given condition to get Fisher's combined test p-value# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ef_Pr_0.001 <- c(pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6)
test.stat_Ef_Pr_0.001 <- -2*sum(log(pvalues_Ef_Pr_0.001))
pfinal_Ef_Pr_0.001<-pchisq(test.stat_Ef_Pr_0.001, df=2*length(pvalues_Ef_Pr_0.001), lower.tail=F)
pvalues_Ef_Pr_0.01 <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6, pvalue_Ef_Pr_0.01_7)
test.stat_Ef_Pr_0.01 <- -2*sum(log(pvalues_Ef_Pr_0.01))
pfinal_Ef_Pr_0.01<-pchisq(test.stat_Ef_Pr_0.01, df=2*length(pvalues_Ef_Pr_0.01), lower.tail=F)
pvalues_Ef_Pr_0.1 <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
test.stat_Ef_Pr_0.1 <- -2*sum(log(pvalues_Ef_Pr_0.1))
pfinal_Ef_Pr_0.1<-pchisq(test.stat_Ef_Pr_0.1, df=2*length(pvalues_Ef_Pr_0.1), lower.tail=F)
pvalues_Ef_Pr_1.0 <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
test.stat_Ef_Pr_1.0 <- -2*sum(log(pvalues_Ef_Pr_1.0))
pfinal_Ef_Pr_1.0<-pchisq(test.stat_Ef_Pr_1.0, df=2*length(pvalues_Ef_Pr_1.0), lower.tail=F)#### Generate Table of final results
pvalues_Ef_Pr_0.001_table <- c(NA,pvalue_Ef_Pr_0.001_2,pvalue_Ef_Pr_0.001_3,pvalue_Ef_Pr_0.001_4,pvalue_Ef_Pr_0.001_5,pvalue_Ef_Pr_0.001_6,NA)
p_Ef_Pr_0.001<-scales::pvalue(pvalues_Ef_Pr_0.001_table)
p_Ef_Pr_0.001 <- p_Ef_Pr_0.001 %>% replace_na("--")
pvalues_Ef_Pr_0.01_table <- c(pvalue_Ef_Pr_0.01_1,pvalue_Ef_Pr_0.01_2,pvalue_Ef_Pr_0.01_3,pvalue_Ef_Pr_0.01_4,pvalue_Ef_Pr_0.01_5,pvalue_Ef_Pr_0.01_6,pvalue_Ef_Pr_0.01_7)
p_Ef_Pr_0.01<-scales::pvalue(pvalues_Ef_Pr_0.01_table)
pvalues_Ef_Pr_0.1_table <- c(pvalue_Ef_Pr_0.1_1,pvalue_Ef_Pr_0.1_2,pvalue_Ef_Pr_0.1_3,pvalue_Ef_Pr_0.1_4,pvalue_Ef_Pr_0.1_5,pvalue_Ef_Pr_0.1_6,pvalue_Ef_Pr_0.1_7)
p_Ef_Pr_0.1<-scales::pvalue(pvalues_Ef_Pr_0.1_table)
pvalues_Ef_Pr_1.0_table <- c(pvalue_Ef_Pr_1.0_1,pvalue_Ef_Pr_1.0_2,pvalue_Ef_Pr_1.0_3,pvalue_Ef_Pr_1.0_4,pvalue_Ef_Pr_1.0_5,pvalue_Ef_Pr_1.0_6,pvalue_Ef_Pr_1.0_7)
p_Ef_Pr_1.0<-scales::pvalue(pvalues_Ef_Pr_1.0_table)
pvalue_table_Ef_Pr <- data.frame("0.001" = p_Ef_Pr_0.001, "0.01" = p_Ef_Pr_0.01, "0.1" = p_Ef_Pr_0.1, "1.0" = p_Ef_Pr_1.0)
colnames(pvalue_table_Ef_Pr) <- c(0.001,0.01, 0.1, 1.0)
pfisher_Ef_Pr<- c(pfinal_Ef_Pr_0.001, pfinal_Ef_Pr_0.01, pfinal_Ef_Pr_0.1, pfinal_Ef_Pr_1.0)
pfisher_Ef_Pr<-scales::pvalue(pfisher_Ef_Pr)
pvalue_table_Ef_Pr[nrow(pvalue_table_Ef_Pr) + 1,]<-pfisher_Ef_Pr
row.names(pvalue_table_Ef_Pr) <- c("Block 1 - June 11th, 2018", "Block 2 - June 22nd, 2018", "Block 3 - June 25th, 2018", "Block 4 - July 19th, 2018", "Block 5 - June 3, 2019", "Block 6 - June 5, 2019", "Block 7 - July 1st, 2019", " Final P-value ")
pvalue_table_Ef_Pr %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
pack_rows("Individual Blocks", 1, 7) %>%
pack_rows("Fisher's Combined", 8, 8)|
Infectious Dose (Abs600nm)
|
||||
|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | |
| Individual Blocks | ||||
| Block 1 - June 11th, 2018 | – | 0.004 | 0.015 | 0.023 |
| Block 2 - June 22nd, 2018 | 0.315 | 0.014 | 0.005 | 0.018 |
| Block 3 - June 25th, 2018 | 0.077 | 0.007 | 0.845 | 0.628 |
| Block 4 - July 19th, 2018 | 0.725 | 0.987 | 0.702 | 0.864 |
| Block 5 - June 3, 2019 | 0.239 | 0.129 | 0.479 | 0.012 |
| Block 6 - June 5, 2019 | 0.010 | 0.184 | 0.620 | 0.872 |
| Block 7 - July 1st, 2019 | – | 0.007 | 0.845 | 0.628 |
| Fisher’s Combined | ||||
| Final P-value | 0.028 | <0.001 | 0.064 | 0.020 |
- Conclusion: Chronic infection with E. faecalis offered significat protection against three out of the four doses. However, the one dose that did not show significant protection (Abs600nm = 0.1) had a borderline p-value = 0.06
Q2 Does having a primary infection influence resistance to secondary infection with P. rettgeri?
Study design: 
- Bacterial load of P. rettgeri was determined at 10 hours post-secondary infection
- For S. marcescens-P. rettgeri, the load was determined by isolating pooled DNA from 5 flies for each condition and using qPCR primer specific for a P. rettgeri gene to assess the relative number of bacteria. This was converted to approximate bacterial load by using a conversion factor based on qPCR and CFU results for flies injected with 20nL of 0.1 Abs600nm of P. rettgeri and then immediately homogenized.
- For E. faecalis - P. rettgeri, the P. rettgeri load was determined using a Colony Forming Unit (CFU) assay using antibiotics to select for P. rettgeri and kill E. faecalis.
- All analyses are with linear modeling followed by anova, with a shapiro test on residuals to assess normality.
- Conclusions
- As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (p<0.0001)
- Chronic infection with S. marcescens and E. faecalis both resulted in lowered calculated bacterial load post secondary infection (p<0.0001).
- There was no significant interaction between chronic infection and secondary infectious dose for either S. marcescens or E. faecalis chronic infection (p=0.30 and 0.78 respectively)
- Residuals pass normality test (Shapiro-Wilk Normality test: p=0.72 and 0.79 respectively)
S. marcescens chronic infection
- The calculated bacterial load after secondary infection with P. rettgeri for flies chronically infected with S. marcescens versus control flies.
- As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (secondary_dose, p<0.0001)
- Primary infection with S. marcescens resulted in lowered calculated bacterial load post secondary infection (primary, p<0.0001).
- There was no significant interaction between primary infection and secondary infectious dose (p=0.30)
- Residuals pass normality test (Shapiro-Wilk Normality test: p=0.72)
Graph
#png(file="Figure3A.png", width=3200, height=1800, res=300)
ggplot(Sm_bacterial_load, aes(x = as.character(secondary_dose), y = log(load_10hr), fill = primary)) +
scale_fill_manual(values=c("white", "#3399FF"),
name=" ",
breaks=c("01_PBS","02_Sm"),
labels=c("control", "chronically infected")) +
geom_boxplot() +
#ggtitle("S. marcescens chronic infection") +
theme_classic(base_size = 24) +
xlab("secondary dose (Abs600nm)") +
ylab("log(calculated bacterial load)") +
guides(fill=guide_legend(title=NULL)) +
geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))
#dev.off()Statistics
S.m_resistance <- lm(log(load_10hr) ~ rep + secondary_dose + primary + secondary_dose:primary, data = Sm_bacterial_load)
S.m_resistance_anova <- anova(S.m_resistance)
S.m_resistance_anova## Analysis of Variance Table
##
## Response: log(load_10hr)
## Df Sum Sq Mean Sq F value Pr(>F)
## rep 8 56.729 7.091 2.8793 0.01059 *
## secondary_dose 1 145.561 145.561 59.1044 6.543e-10 ***
## primary 1 94.077 94.077 38.1995 1.328e-07 ***
## secondary_dose:primary 1 2.651 2.651 1.0765 0.30467
## Residuals 48 118.213 2.463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(resid(S.m_resistance))##
## Shapiro-Wilk normality test
##
## data: resid(S.m_resistance)
## W = 0.98595, p-value = 0.7192E. faecalis chronic infection
- The bacterial load after secondary infection with P. rettgeri for flies chronically infected with E. faecalis versus control flies.
- As expected based on previously published results (Chambers et al 2019), infectious dose of the secondary P. rettgeri infection significantly impacted bacterial load 10 hours after infection (secondary_dose, p<0.0001)
- Primary infection with E. faecalis resulted in lowered CFU/fly post secondary infection (primary, p = 0.0001).
- There was no significant interaction between primary infection and secondary infectious dose (p=0.78)
- Residuals pass normality test (Shapiro-Wilk Normality test: p=0.79)
Graph
#BACTERIAL LOAD
Ef_bacterial_load_10hr <- subset(Ef_bacterial_load, time == 10)
#png(file="Figure3B.png",width=3200, height=1800, res=300)
ggplot(Ef_bacterial_load_10hr, aes(x = as.character(secondary_dose), y = log(load), fill = primary)) +
scale_fill_manual(values=c("white", "#3399FF"),
name=" ",
breaks=c("01_PBS","02_Efae"),
labels=c("control", "chronically infected")) +
geom_boxplot() +
#ggtitle("E. faecalis chronic infection") +
theme_classic(base_size = 24) +
xlab("secondary dose (Abs600nm)") +
ylab("log(CFU/fly)") +
guides(fill=guide_legend(title=NULL)) +
geom_point(color = "black", size = 1.3, position=position_jitterdodge(jitter.width=0.2))
#dev.off()Statistics
E.fae_resistance <- lm(log(load) ~ replicate + secondary_dose + primary + secondary_dose:primary, data = Ef_bacterial_load_10hr)
E.fae_resistance_anova <- anova(E.fae_resistance)
E.fae_resistance_anova## Analysis of Variance Table
##
## Response: log(load)
## Df Sum Sq Mean Sq F value Pr(>F)
## replicate 5 231.74 46.35 15.403 3.504e-13 ***
## secondary_dose 1 799.17 799.17 265.589 < 2.2e-16 ***
## primary 1 35.09 35.09 11.660 0.0007463 ***
## secondary_dose:primary 1 0.23 0.23 0.075 0.7843830
## Residuals 247 743.23 3.01
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shapiro.test(resid(E.fae_resistance))##
## Shapiro-Wilk normality test
##
## data: resid(E.fae_resistance)
## W = 0.99623, p-value = 0.7984Figure 3
k <- ggdraw() +
ggtitle('P. rettgeri') +
draw_image("Figure3A.png")
l <- ggdraw() +
draw_image("Figure3B.png")
m <- ggdraw() +
draw_image("Figure3C.png")
n <- ggdraw() +
draw_image("Figure3D.png")
multiplot3 <- k + m + l + n
x<-multiplot3
#plot_annotation(tag_levels = 'A')
ggsave("Fig.3.ResistanceTolerance.pdf", x, device = "pdf")
x
Q3 Does the protection conferred by primary infection with S. marcescens protect against infections that cause very high mortality?
Study design: 
- Since we saw such robust protection with S. marcescens, we wanted to determine whether it also protected against a more lethal infection with P. sneebia.
- Initially analyzed using Cox Proportional Hazards mixed effects
model with date of injection, primary injection and secondary dose as
fixed effects and random effect terms accounting for experimental
structure (e.g. date of injection, housing of groups of flies in vials).
However, this model violated the assumptions of Cox Proportional
Hazards. + Therefore we switched to using log-rank pairwise comparisons
to determine whether chronic infection significantly impacted survival
of secondary infection at an individual dose.
- LogRank comparisons were done on data from each date individually
- If experiments on all dates showed higher survival in conditions with chronic infection, p-values were combined using Fisher’s test and only the final p-value is reported in the text.
- If there was a variable effect of chronic infection, the number of
individual experiments with significant effects is reported.
- Conclusions
- S. marcescens chronic infection provided significant protection during secondary infection with P. sneebia at all infectious doses (p<0.001)
Total Graph
#table(Sm_Ps$primary, Sm_Ps$secondary_dose)
#png(file="Sm_Ps.png", width = 1700, height = 1500, res = 300)
#par(mar=c(4,4,0.1,0.1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Sm_Ps$hour, Sm_Ps$status) ~ Sm_Ps$primary + Sm_Ps$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.6, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "PBS-0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Appendix Figure 1B - Controls
Sm_Ps_PBS<-subset(Sm_Ps,secondary_dose=="0")
#table(Sm_Ps_PBS$primary)
#THE GRAPH
#png(file="AFigure1B.png",width=1200,height=900,res=130)
par(mar=c(4,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(day, status) ~ primary + secondary_dose, data=Sm_Ps_PBS), col=c("gray"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("double injection control (n=59)", "chronic infection only (n=59)"),col=c("gray"),lty=c(3, 1),lwd=5, cex=1)
#dev.off()Individual Figure Panels
Figure 4A - 0.001
Sm_Ps_0.001<-subset(Sm_Ps,secondary_dose=="0.001")
#table(Sm_Ps_0.001$primary)
#THE GRAPH
#png(file="Figure4A.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.001), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=82)", "chronic & secondary infection (n=66)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()Figure 4B - 0.01
Sm_Ps_0.01<-subset(Sm_Ps,secondary_dose=="0.01")
#table(Sm_Ps_0.01$primary)
#THE GRAPH
#png(file="Figure4B.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.01), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("bottomleft",inset=0, bty="n" , c("secondary infection only (n=87)", "chronic & secondary infection (n=80)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()Figure 4C - 0.1
Sm_Ps_0.1<-subset(Sm_Ps,secondary_dose=="0.1")
#table(Sm_Ps_0.1$primary)
#THE GRAPH
#png(file="Figure4C.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_0.1), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=80)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()Figure 4D - 1.0
Sm_Ps_1.0<-subset(Sm_Ps,secondary_dose=="1")
#table(Sm_Ps_1.0$primary)
#png(file="Figure4D.png",width=1200,height=900,res=130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(hour, status) ~ primary + secondary_dose, data=Sm_Ps_1.0), col=c("#003694"), lty=c(3, 1), lwd=5,yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2, xlim=c(0,72), xaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
#THE LEGEND
legend("topright",inset=0, bty="n" , c("secondary infection only (n=75)", "chronic & secondary infection (n=81)"),col=c("#003694"),lty=c(3, 1),lwd=5, cex=1.6)
#dev.off()Figure 4
m <- ggdraw() +
draw_image("Figure4A.png")
n <- ggdraw() +
draw_image("Figure4B.png")
o <- ggdraw() +
ggtitle('P. rettgeri') +
draw_image("Figure4C.png")
p <- ggdraw() +
draw_image("Figure4D.png")
multiplot3 <- m + n + o + p
w<-multiplot3 +
plot_annotation(tag_levels = 'A')
ggsave("Fig.4.Psneebia.pdf", w, device = "pdf")
w
Individual Blocks
Block 1
Sm_Ps_block1 <- subset(Sm_Ps, date == "10/1/2020")
table(Sm_Ps_block1$primary, Sm_Ps_block1$secondary_dose)##
## 0 0.001 0.01 0.1 1
## PBS 19 30 29 20 30
## S.m 19 16 24 29 25#subset by secondary dose
Sm_Ps_0.001_1<-subset(Sm_Ps_block1, secondary_dose=="0.001")
Sm_Ps_0.01_1<-subset(Sm_Ps_block1, secondary_dose=="0.01")
Sm_Ps_0.1_1<-subset(Sm_Ps_block1, secondary_dose=="0.1")
Sm_Ps_1.0_1<-subset(Sm_Ps_block1, secondary_dose=="1")
surv_Sm_Ps_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_1)
surv_Sm_Ps_0.001_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 21 18.3 0.388 1.13
## primary=S.m 16 10 12.7 0.562 1.13
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3surv_Sm_Ps_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_1)
surv_Sm_Ps_0.01_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 29 29 14.7 13.87 29.7
## primary=S.m 24 23 37.3 5.47 29.7
##
## Chisq= 29.7 on 1 degrees of freedom, p= 5e-08surv_Sm_Ps_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_1)
surv_Sm_Ps_0.1_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 20 20 6.8 25.65 45.2
## primary=S.m 29 29 42.2 4.13 45.2
##
## Chisq= 45.2 on 1 degrees of freedom, p= 2e-11surv_Sm_Ps_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_1)
surv_Sm_Ps_1.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 18.3 7.51 19.9
## primary=S.m 25 25 36.7 3.74 19.9
##
## Chisq= 19.9 on 1 degrees of freedom, p= 8e-06plot(survfit(Surv(Sm_Ps_block1$hour, Sm_Ps_block1$status) ~ Sm_Ps_block1$primary + Sm_Ps_block1$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Block 2
Sm_Ps_block2 <- subset(Sm_Ps, date == "10/10/2020")
table(Sm_Ps_block2$primary, Sm_Ps_block2$secondary_dose)##
## 0 0.001 0.01 0.1 1
## PBS 20 30 30 30 30
## S.m 20 30 30 30 30#subset by secondary dose
Sm_Ps_0.001_2<-subset(Sm_Ps_block2, secondary_dose=="0.001")
Sm_Ps_0.01_2<-subset(Sm_Ps_block2, secondary_dose=="0.01")
Sm_Ps_0.1_2<-subset(Sm_Ps_block2, secondary_dose=="0.1")
Sm_Ps_1.0_2<-subset(Sm_Ps_block2, secondary_dose=="1")
surv_Sm_Ps_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_2)
surv_Sm_Ps_0.001_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 25 16 5.07 10.5
## primary=S.m 30 13 22 3.69 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001surv_Sm_Ps_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_2)
surv_Sm_Ps_0.01_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 15.8 12.77 24
## primary=S.m 30 26 40.2 5.02 24
##
## Chisq= 24 on 1 degrees of freedom, p= 9e-07surv_Sm_Ps_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_2)
surv_Sm_Ps_0.1_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 14.6 16.35 40.9
## primary=S.m 30 30 45.4 5.24 40.9
##
## Chisq= 40.9 on 1 degrees of freedom, p= 2e-10surv_Sm_Ps_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_2)
surv_Sm_Ps_1.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 13.5 20.14 44
## primary=S.m 30 30 46.5 5.85 44
##
## Chisq= 44 on 1 degrees of freedom, p= 3e-11plot(survfit(Surv(Sm_Ps_block2$hour, Sm_Ps_block2$status) ~ Sm_Ps_block2$primary + Sm_Ps_block2$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Block 3
Sm_Ps_block3 <- subset(Sm_Ps, date == "10/24/2020")
table(Sm_Ps_block3$primary, Sm_Ps_block3$secondary_dose)##
## 0 0.001 0.01 0.1 1
## PBS 20 22 28 30 15
## S.m 19 20 26 22 26#subset by secondary dose
Sm_Ps_0.001_3<-subset(Sm_Ps_block3, secondary_dose=="0.001")
Sm_Ps_0.01_3<-subset(Sm_Ps_block3, secondary_dose=="0.01")
Sm_Ps_0.1_3<-subset(Sm_Ps_block3, secondary_dose=="0.1")
Sm_Ps_1.0_3<-subset(Sm_Ps_block3, secondary_dose=="1")
surv_Sm_Ps_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.001_3)
surv_Sm_Ps_0.001_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 22 21 15.6 1.90 9.75
## primary=S.m 20 13 18.4 1.61 9.75
##
## Chisq= 9.8 on 1 degrees of freedom, p= 0.002surv_Sm_Ps_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.01_3)
surv_Sm_Ps_0.01_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 28 28 13.1 16.80 34.7
## primary=S.m 26 21 35.9 6.16 34.7
##
## Chisq= 34.7 on 1 degrees of freedom, p= 4e-09surv_Sm_Ps_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_0.1_3)
surv_Sm_Ps_0.1_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 30 30 14 18.33 40.2
## primary=S.m 22 21 37 6.93 40.2
##
## Chisq= 40.2 on 1 degrees of freedom, p= 2e-10surv_Sm_Ps_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Sm_Ps_1.0_3)
surv_Sm_Ps_1.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Sm_Ps_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 15 15 4.8 21.71 34.2
## primary=S.m 26 26 36.2 2.88 34.2
##
## Chisq= 34.2 on 1 degrees of freedom, p= 5e-09plot(survfit(Surv(Sm_Ps_block3$hour, Sm_Ps_block3$status) ~ Sm_Ps_block3$primary + Sm_Ps_block3$secondary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1), lwd=2.3, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-PBS", "0.001", "PBS-0.01", "PBS-0.1", "PBS-1.0", "S.m-PBS", "S.m-0.001", "S.m-0.01", "S.m-0.1", "S.m-1.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5"), lty=c(3, 3, 3, 3, 3, 1, 1, 1, 1, 1),lwd=4, cex=1)
Statistical Tests
- In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Sm_Ps_0.001_1<-pchisq(surv_Sm_Ps_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_1<-pchisq(surv_Sm_Ps_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_1<-pchisq(surv_Sm_Ps_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_1<-pchisq(surv_Sm_Ps_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.001_2<-pchisq(surv_Sm_Ps_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_2<-pchisq(surv_Sm_Ps_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_2<-pchisq(surv_Sm_Ps_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_2<-pchisq(surv_Sm_Ps_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.001_3<-pchisq(surv_Sm_Ps_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.01_3<-pchisq(surv_Sm_Ps_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_0.1_3<-pchisq(surv_Sm_Ps_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Sm_Ps_1.0_3<-pchisq(surv_Sm_Ps_1.0_3$chisq, df=1, lower.tail=F)# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Sm_Ps_0.001 <- c(pvalue_Sm_Ps_0.001_2,pvalue_Sm_Ps_0.001_3)
test.stat_Sm_Ps_0.001 <- -2*sum(log(pvalues_Sm_Ps_0.001))
pfinal_Sm_Ps_0.001<-pchisq(test.stat_Sm_Ps_0.001, df=2*length(pvalues_Sm_Ps_0.001), lower.tail=F)
pvalues_Sm_Ps_0.01 <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
test.stat_Sm_Ps_0.01 <- -2*sum(log(pvalues_Sm_Ps_0.01))
pfinal_Sm_Ps_0.01<-pchisq(test.stat_Sm_Ps_0.01, df=2*length(pvalues_Sm_Ps_0.01), lower.tail=F)
pvalues_Sm_Ps_0.1 <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
test.stat_Sm_Ps_0.1 <- -2*sum(log(pvalues_Sm_Ps_0.1))
pfinal_Sm_Ps_0.1<-pchisq(test.stat_Sm_Ps_0.1, df=2*length(pvalues_Sm_Ps_0.1), lower.tail=F)
pvalues_Sm_Ps_1.0 <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
test.stat_Sm_Ps_1.0 <- -2*sum(log(pvalues_Sm_Ps_1.0))
pfinal_Sm_Ps_1.0<-pchisq(test.stat_Sm_Ps_1.0, df=2*length(pvalues_Sm_Ps_1.0), lower.tail=F)pvalues_Sm_Ps_0.001_table <- c(pvalue_Sm_Ps_0.001_1, pvalue_Sm_Ps_0.001_2, pvalue_Sm_Ps_0.001_3)
p_Sm_Ps_0.001<-scales::pvalue(pvalues_Sm_Ps_0.001_table)
pvalues_Sm_Ps_0.01_table <- c(pvalue_Sm_Ps_0.01_1,pvalue_Sm_Ps_0.01_2,pvalue_Sm_Ps_0.01_3)
p_Sm_Ps_0.01<-scales::pvalue(pvalues_Sm_Ps_0.01_table)
pvalues_Sm_Ps_0.1_table <- c(pvalue_Sm_Ps_0.1_1,pvalue_Sm_Ps_0.1_2,pvalue_Sm_Ps_0.1_3)
p_Sm_Ps_0.1<-scales::pvalue(pvalues_Sm_Ps_0.1_table)
pvalues_Sm_Ps_1.0_table <- c(pvalue_Sm_Ps_1.0_1,pvalue_Sm_Ps_1.0_2,pvalue_Sm_Ps_1.0_3)
p_Sm_Ps_1.0<-scales::pvalue(pvalues_Sm_Ps_1.0_table)
pvalue_table_Sm_Ps <- data.frame("0.001" = p_Sm_Ps_0.001, "0.01" = p_Sm_Ps_0.01, "0.1" = p_Sm_Ps_0.1, "1.0" = p_Sm_Ps_1.0)
colnames(pvalue_table_Sm_Ps) <- c(0.001,0.01, 0.1, 1.0)
pfisher_Sm_Ps<- c(pfinal_Sm_Ps_0.001, pfinal_Sm_Ps_0.01, pfinal_Sm_Ps_0.1, pfinal_Sm_Ps_1.0)
pfisher_Sm_Ps<-scales::pvalue(pfisher_Sm_Ps)
pvalue_table_Sm_Ps[nrow(pvalue_table_Sm_Ps) + 1,]<-pfisher_Sm_Ps
row.names(pvalue_table_Sm_Ps) <- c("Block 1 - October 1st, 2021", "Block 2 - October 10th, 2021", "Block 3 - October 24th, 2021", " Final P-value ")
pvalue_table_Sm_Ps %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Infectious Dose (Abs600nm)" = 4)) %>%
pack_rows("Individual Blocks", 1, 3) %>%
pack_rows("Fisher's Combined", 4, 4)|
Infectious Dose (Abs600nm)
|
||||
|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | |
| Individual Blocks | ||||
| Block 1 - October 1st, 2021 | 0.287 | <0.001 | <0.001 | <0.001 |
| Block 2 - October 10th, 2021 | 0.001 | <0.001 | <0.001 | <0.001 |
| Block 3 - October 24th, 2021 | 0.002 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | ||||
| Final P-value | <0.001 | <0.001 | <0.001 | <0.001 |
- Conclusion: Chronic infection with S. marcescens offered significat protection against all doses of P. sneebia (p<0.001)
Q4 Does the dose of the primary infection impact the strength of protection during secondary infection?
- The infectious dose impacts the number of bacteria present in the fly during chronic infection (Chambers et al. 2019), therefore we hypothesized that this might influence the strength of protection during secondary infection.
- Flies were injected with one of 5 doses of S. marcescens during primary infection (0.001, 0.01, 0.1, 1.0 and 2.0) or PBS. Then they were given a challenge dose of either P. rettgeri or P. sneebia (2.0 Abs600nM) as a secondary infection.
Study design: 
+ Initially
analyzed using Cox Proportional Hazards mixed effects model with date of
injection, primary injection dose as fixed effects and random effect
terms accounting for experimental structure (e.g. date of injection,
housing of groups of flies in vials). However, this model violated the
assumptions of Cox Proportional Hazards. + Therefore we switched to
using log-rank pairwise comparisons to determine whether dose used to
induce chronic infection significantly impacted survival of secondary
infection by comparing to the control flies given only a secondary
infection. + LogRank comparisons were done on data from each date
individually + If experiments on all dates showed higher survival in
conditions with chronic infection, p-values were combined using Fisher’s
test and only the final p-value is reported in the text. + If there was
a variable effect of chronic infection, the number of individual
experiments with significant effects is reported. +
Conclusions + Moderate doses of S. marcescens
(Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection
against secondary infection with both P. rettgeri ad P.
sneebia (P<0.001) + The lowest dose of S. marcescens
(Abs600nm = 0.01, 0.1 & 1.0) did not provide protection against
either secondary infection. + The highest dose provided inconsistent
protection, sometimes providing robust protection and other times dying
as quickly as (or even faster than) control flies carrying only a
secondary infection.
P. rettgeri secondary infection
- Conclusion: Moderate doses of the primary infection (Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection against infection (P<0.001).
- The lowest dose (Abs600nm = 0.001) did not provide significant protection
- The highest dose provided inconsistent protection with 8 out of the 9 blocks demonstrating better survival than controls and one with worse survival than controls. In Block 1, flies given a primary infection at the highest dose died more than controls with only the secondary infection.
Figure 5A
Pr_varied_primary = subset(Pr_varied_primary, secondary == "P.rett")
Pr_varied_primary$primary_dose <- as.factor(Pr_varied_primary$primary_dose)
table(Pr_varied_primary$date, Pr_varied_primary$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 10/18/2020 40 38 35 40 27 31
## 10/20/2020 36 37 35 39 35 30
## 10/29/2020 40 37 46 36 37 33
## 11/2/2020 37 31 30 40 30 27
## 11/28/2020 39 47 47 48 45 42
## 11/7/2020 40 38 40 35 34 36
## 12/4/2020 35 38 36 34 27 0
## 6/17/2021 37 47 47 44 42 36
## 6/18/2021 33 48 48 42 42 63table(Pr_varied_primary$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 337 361 364 358 319 298#png(file="Figure5A.png", width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Pr_varied_primary$day, Pr_varied_primary$status) ~ Pr_varied_primary$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="days post-secondary infection", cex.lab=2.5, xlim=c(0,7), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=4.5)
box(lwd=5)
#legend("bottomleft", bty="n" , c("secondary infection only", "chronic infection (0.001)", "chronic infection (0.01)", "chronic infection (0.1)", "chronic infection (1.0)", "chronic infection (2.0)"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)
#dev.off()Individual blocks
Block 1
Pr_varied_primary_block1 <- subset(Pr_varied_primary, date == "10/18/2020")
table(Pr_varied_primary_block1$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 40 38 35 40 27 31#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_1<-subset(Pr_varied_primary_block1, primary_dose == "0.001"|primary_dose=="0")
table(Pr_varied_primary_0.001_1$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 40 38 0 0 0 0Pr_varied_primary_0.01_1<-subset(Pr_varied_primary_block1, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_1<-subset(Pr_varied_primary_block1, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_1<-subset(Pr_varied_primary_block1, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_1<-subset(Pr_varied_primary_block1, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_1)
surv_Pr_varied_primary_0.001_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 37.4 0.792 3.24
## primary=S.m 38 35 29.6 1.003 3.24
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07surv_Pr_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_1)
surv_Pr_varied_primary_0.01_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 28.8 0.358 1.11
## primary=S.m 35 24 27.2 0.379 1.11
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3surv_Pr_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_1)
surv_Pr_varied_primary_0.1_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 18.6 9.60 22.5
## primary=S.m 40 12 25.4 7.05 22.5
##
## Chisq= 22.5 on 1 degrees of freedom, p= 2e-06surv_Pr_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_1)
surv_Pr_varied_primary_1.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 26.3 1.22 3.27
## primary=S.m 27 22 27.7 1.16 3.27
##
## Chisq= 3.3 on 1 degrees of freedom, p= 0.07surv_Pr_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_1)
surv_Pr_varied_primary_2.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 32 37.4 0.776 4.16
## primary=S.m 31 31 25.6 1.133 4.16
##
## Chisq= 4.2 on 1 degrees of freedom, p= 0.04par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block1$day, Pr_varied_primary_block1$status) ~ Pr_varied_primary_block1$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 2
Pr_varied_primary_block2 <- subset(Pr_varied_primary, date == "10/20/2020")
table(Pr_varied_primary_block2$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 36 37 35 39 35 30#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_2<-subset(Pr_varied_primary_block2, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_2<-subset(Pr_varied_primary_block2, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_2<-subset(Pr_varied_primary_block2, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_2<-subset(Pr_varied_primary_block2, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_2<-subset(Pr_varied_primary_block2, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_2)
surv_Pr_varied_primary_0.001_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 32.6 0.0783 0.403
## primary=S.m 37 35 33.4 0.0764 0.403
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5surv_Pr_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_2)
surv_Pr_varied_primary_0.01_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 27.7 0.388 1.57
## primary=S.m 35 28 31.3 0.344 1.57
##
## Chisq= 1.6 on 1 degrees of freedom, p= 0.2surv_Pr_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_2)
surv_Pr_varied_primary_0.1_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 19.6 6.63 16.7
## primary=S.m 39 18 29.4 4.42 16.7
##
## Chisq= 16.7 on 1 degrees of freedom, p= 4e-05surv_Pr_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_2)
surv_Pr_varied_primary_1.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 21.7 4.03 9.17
## primary=S.m 35 28 37.3 2.34 9.17
##
## Chisq= 9.2 on 1 degrees of freedom, p= 0.002surv_Pr_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_2)
surv_Pr_varied_primary_2.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 36 31 25.5 1.192 3.23
## primary=S.m 30 27 32.5 0.935 3.23
##
## Chisq= 3.2 on 1 degrees of freedom, p= 0.07par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block2$day, Pr_varied_primary_block2$status) ~ Pr_varied_primary_block2$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 3
Pr_varied_primary_block3 <- subset(Pr_varied_primary, date == "10/29/2020")
table(Pr_varied_primary_block3$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 40 37 46 36 37 33#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_3<-subset(Pr_varied_primary_block3, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_3<-subset(Pr_varied_primary_block3, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_3<-subset(Pr_varied_primary_block3, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_3<-subset(Pr_varied_primary_block3, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_3<-subset(Pr_varied_primary_block3, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_3)
surv_Pr_varied_primary_0.001_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 34.6 0.0559 0.235
## primary=S.m 37 31 32.4 0.0597 0.235
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6surv_Pr_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_3)
surv_Pr_varied_primary_0.01_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 28.6 1.93 5.61
## primary=S.m 46 32 39.4 1.40 5.61
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02surv_Pr_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_3)
surv_Pr_varied_primary_0.1_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 18.3 17.1 41.4
## primary=S.m 36 9 26.7 11.7 41.4
##
## Chisq= 41.4 on 1 degrees of freedom, p= 1e-10surv_Pr_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_3)
surv_Pr_varied_primary_1.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 19.6 13.83 30.3
## primary=S.m 37 19 35.4 7.63 30.3
##
## Chisq= 30.3 on 1 degrees of freedom, p= 4e-08surv_Pr_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_3)
surv_Pr_varied_primary_2.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 36 27.8 2.43 6.33
## primary=S.m 33 27 35.2 1.92 6.33
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block3$day, Pr_varied_primary_block3$status) ~ Pr_varied_primary_block3$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 4
Pr_varied_primary_block4 <- subset(Pr_varied_primary, date == "11/2/2020")
table(Pr_varied_primary_block4$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 37 31 30 40 30 27#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_4<-subset(Pr_varied_primary_block4, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_4<-subset(Pr_varied_primary_block4, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_4<-subset(Pr_varied_primary_block4, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_4<-subset(Pr_varied_primary_block4, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_4<-subset(Pr_varied_primary_block4, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_4)
surv_Pr_varied_primary_0.001_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 33.2 0.000817 0.00537
## primary=S.m 31 28 27.8 0.000973 0.00537
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9surv_Pr_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_4)
surv_Pr_varied_primary_0.01_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 26.5 1.57 6.91
## primary=S.m 30 20 26.5 1.58 6.91
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.009surv_Pr_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_4)
surv_Pr_varied_primary_0.1_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 21.2 6.58 20.5
## primary=S.m 40 19 30.8 4.52 20.5
##
## Chisq= 20.5 on 1 degrees of freedom, p= 6e-06surv_Pr_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_4)
surv_Pr_varied_primary_1.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 20.9 7.08 19.2
## primary=S.m 30 19 31.1 4.74 19.2
##
## Chisq= 19.2 on 1 degrees of freedom, p= 1e-05surv_Pr_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_4)
surv_Pr_varied_primary_2.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 33 30.6 0.181 0.686
## primary=S.m 27 27 29.4 0.189 0.686
##
## Chisq= 0.7 on 1 degrees of freedom, p= 0.4par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block4$day, Pr_varied_primary_block4$status) ~ Pr_varied_primary_block4$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 5
Pr_varied_primary_block5 <- subset(Pr_varied_primary, date == "11/7/2020")
table(Pr_varied_primary_block5$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 40 38 40 35 34 36#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_5<-subset(Pr_varied_primary_block5, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_5<-subset(Pr_varied_primary_block5, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_5<-subset(Pr_varied_primary_block5, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_5<-subset(Pr_varied_primary_block5, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_5<-subset(Pr_varied_primary_block5, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_5)
surv_Pr_varied_primary_0.001_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 33.9 0.0383 0.156
## primary=S.m 38 32 33.1 0.0391 0.156
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7surv_Pr_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_5)
surv_Pr_varied_primary_0.01_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 32 0.289 0.942
## primary=S.m 40 35 38 0.243 0.942
##
## Chisq= 0.9 on 1 degrees of freedom, p= 0.3surv_Pr_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_5)
surv_Pr_varied_primary_0.1_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 21.5 8.41 21.9
## primary=S.m 35 15 28.5 6.37 21.9
##
## Chisq= 21.9 on 1 degrees of freedom, p= 3e-06surv_Pr_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_5)
surv_Pr_varied_primary_1.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 20 11.28 27.4
## primary=S.m 34 16 31 7.27 27.4
##
## Chisq= 27.4 on 1 degrees of freedom, p= 2e-07surv_Pr_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_5)
surv_Pr_varied_primary_2.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 40 35 26.1 3.02 6.7
## primary=S.m 36 34 42.9 1.84 6.7
##
## Chisq= 6.7 on 1 degrees of freedom, p= 0.01par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block5$day, Pr_varied_primary_block5$status) ~ Pr_varied_primary_block5$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 6
Pr_varied_primary_block6 <- subset(Pr_varied_primary, date == "11/28/2020")
table(Pr_varied_primary_block6$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 39 47 47 48 45 42#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_6<-subset(Pr_varied_primary_block6, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_6<-subset(Pr_varied_primary_block6, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_6<-subset(Pr_varied_primary_block6, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_6<-subset(Pr_varied_primary_block6, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_6<-subset(Pr_varied_primary_block6, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_6)
surv_Pr_varied_primary_0.001_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 38.1 0.679 3.44
## primary=S.m 47 44 38.9 0.664 3.44
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06surv_Pr_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_6)
surv_Pr_varied_primary_0.01_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 32.7 0.00346 0.0119
## primary=S.m 47 39 39.3 0.00288 0.0119
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9surv_Pr_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_6)
surv_Pr_varied_primary_0.1_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 19.6 9.14 19.7
## primary=S.m 48 23 36.4 4.93 19.7
##
## Chisq= 19.7 on 1 degrees of freedom, p= 9e-06surv_Pr_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_6)
surv_Pr_varied_primary_1.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 19.1 10.03 21.9
## primary=S.m 45 20 33.9 5.67 21.9
##
## Chisq= 21.9 on 1 degrees of freedom, p= 3e-06surv_Pr_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_6)
surv_Pr_varied_primary_2.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 33 25.8 2.03 4.56
## primary=S.m 42 37 44.2 1.18 4.56
##
## Chisq= 4.6 on 1 degrees of freedom, p= 0.03par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block6$day, Pr_varied_primary_block6$status) ~ Pr_varied_primary_block6$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 7
Pr_varied_primary_block7 <- subset(Pr_varied_primary, date == "12/4/2020")
table(Pr_varied_primary_block7$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 35 38 36 34 27 0#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_7<-subset(Pr_varied_primary_block7, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_7<-subset(Pr_varied_primary_block7, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_7<-subset(Pr_varied_primary_block7, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_7<-subset(Pr_varied_primary_block7, primary_dose=="1"|primary_dose=="0")
#Pr_varied_primary_2.0_7<-subset(Pr_varied_primary_block7, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_7)
surv_Pr_varied_primary_0.001_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 31.1 0.000230 0.00162
## primary=S.m 38 33 32.9 0.000217 0.00162
##
## Chisq= 0 on 1 degrees of freedom, p= 1surv_Pr_varied_primary_0.01_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_7)
surv_Pr_varied_primary_0.01_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 24.8 1.53 7.31
## primary=S.m 36 21 27.2 1.40 7.31
##
## Chisq= 7.3 on 1 degrees of freedom, p= 0.007surv_Pr_varied_primary_0.1_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_7)
surv_Pr_varied_primary_0.1_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 16.1 13.8 40.8
## primary=S.m 34 6 20.9 10.6 40.8
##
## Chisq= 40.8 on 1 degrees of freedom, p= 2e-10surv_Pr_varied_primary_1.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_7)
surv_Pr_varied_primary_1.0_7## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_7)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 31 20.8 4.96 14.1
## primary=S.m 27 21 31.2 3.31 14.1
##
## Chisq= 14.1 on 1 degrees of freedom, p= 2e-04#surv_Pr_varied_primary_2.0_7<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_7)
#surv_Pr_varied_primary_2.0_7
par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block7$day, Pr_varied_primary_block7$status) ~ Pr_varied_primary_block7$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 8
Pr_varied_primary_block8 <- subset(Pr_varied_primary, date == "6/17/2021")
table(Pr_varied_primary_block8$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 37 47 47 44 42 36#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_8<-subset(Pr_varied_primary_block8, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_8<-subset(Pr_varied_primary_block8, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_8<-subset(Pr_varied_primary_block8, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_8<-subset(Pr_varied_primary_block8, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_8<-subset(Pr_varied_primary_block8, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_8)
surv_Pr_varied_primary_0.001_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 36.2 0.0185 0.123
## primary=S.m 47 44 44.8 0.0150 0.123
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7surv_Pr_varied_primary_0.01_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_8)
surv_Pr_varied_primary_0.01_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 32.2 0.700 3.45
## primary=S.m 47 42 46.8 0.483 3.45
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.06surv_Pr_varied_primary_0.1_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_8)
surv_Pr_varied_primary_0.1_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 17.5 21.9 55.4
## primary=S.m 44 14 33.5 11.4 55.4
##
## Chisq= 55.4 on 1 degrees of freedom, p= 1e-13surv_Pr_varied_primary_1.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_8)
surv_Pr_varied_primary_1.0_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 16 27.4 67.1
## primary=S.m 42 13 34 12.9 67.1
##
## Chisq= 67.1 on 1 degrees of freedom, p= 3e-16surv_Pr_varied_primary_2.0_8<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_8)
surv_Pr_varied_primary_2.0_8## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_8)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 24 7.02 27.6
## primary=S.m 36 35 48 3.51 27.6
##
## Chisq= 27.6 on 1 degrees of freedom, p= 2e-07par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block8$day, Pr_varied_primary_block8$status) ~ Pr_varied_primary_block8$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 9
Pr_varied_primary_block9<- subset(Pr_varied_primary, date == "6/18/2021")
table(Pr_varied_primary_block9$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 33 48 48 42 42 63#subset by primary dose for pairwise comparison dose
Pr_varied_primary_0.001_9<-subset(Pr_varied_primary_block9, primary_dose == "0.001"|primary_dose=="0")
Pr_varied_primary_0.01_9<-subset(Pr_varied_primary_block9, primary_dose=="0.01"|primary_dose=="0")
Pr_varied_primary_0.1_9<-subset(Pr_varied_primary_block9, primary_dose=="0.1"|primary_dose=="0")
Pr_varied_primary_1.0_9<-subset(Pr_varied_primary_block9, primary_dose=="1"|primary_dose=="0")
Pr_varied_primary_2.0_9<-subset(Pr_varied_primary_block9, primary_dose=="2"|primary_dose=="0")
surv_Pr_varied_primary_0.001_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.001_9)
surv_Pr_varied_primary_0.001_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.001_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 32.2 0.0206 1.39
## primary=S.m 48 46 46.8 0.0142 1.39
##
## Chisq= 1.4 on 1 degrees of freedom, p= 0.2surv_Pr_varied_primary_0.01_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.01_9)
surv_Pr_varied_primary_0.01_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.01_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 27.7 1.013 10.5
## primary=S.m 48 38 43.3 0.648 10.5
##
## Chisq= 10.5 on 1 degrees of freedom, p= 0.001surv_Pr_varied_primary_0.1_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_0.1_9)
surv_Pr_varied_primary_0.1_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_0.1_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 19.8 8.80 38.8
## primary=S.m 42 18 31.2 5.58 38.8
##
## Chisq= 38.8 on 1 degrees of freedom, p= 5e-10surv_Pr_varied_primary_1.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_1.0_9)
surv_Pr_varied_primary_1.0_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_1.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 17.2 14.6 53.7
## primary=S.m 42 13 28.8 8.7 53.7
##
## Chisq= 53.7 on 1 degrees of freedom, p= 2e-13surv_Pr_varied_primary_2.0_9<-survdiff(Surv(day, status) ~ primary, data=Pr_varied_primary_2.0_9)
surv_Pr_varied_primary_2.0_9## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Pr_varied_primary_2.0_9)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 33 33 19.2 9.82 35.5
## primary=S.m 63 50 63.8 2.97 35.5
##
## Chisq= 35.5 on 1 degrees of freedom, p= 2e-09par(mar=c(4,4,0.1,0.1))
plot(survfit(Surv(Pr_varied_primary_block9$day, Pr_varied_primary_block9$status) ~ Pr_varied_primary_block9$primary_dose), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Days Post-Secondary Infection", cex.lab=1.3, xlim=c(0,7), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.rett 2.0", "S.m 0.001-P.rett 2.0", "S.m 0.01-P.rett 2.0", "S.m 0.1-P.rett 2.0", "S.m 1.0-P.rett 2.0", "S.m 2.0-P.rett 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Statistical Tests
- In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Pr_varied_primary_0.001_1<-pchisq(surv_Pr_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_1<-pchisq(surv_Pr_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_1<-pchisq(surv_Pr_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_1<-pchisq(surv_Pr_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_1<-pchisq(surv_Pr_varied_primary_2.0_1$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_2<-pchisq(surv_Pr_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_2<-pchisq(surv_Pr_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_2<-pchisq(surv_Pr_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_2<-pchisq(surv_Pr_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_2<-pchisq(surv_Pr_varied_primary_2.0_2$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_3<-pchisq(surv_Pr_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_3<-pchisq(surv_Pr_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_3<-pchisq(surv_Pr_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_3<-pchisq(surv_Pr_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_3<-pchisq(surv_Pr_varied_primary_2.0_3$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_4<-pchisq(surv_Pr_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_4<-pchisq(surv_Pr_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_4<-pchisq(surv_Pr_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_4<-pchisq(surv_Pr_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_4<-pchisq(surv_Pr_varied_primary_2.0_4$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_5<-pchisq(surv_Pr_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_5<-pchisq(surv_Pr_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_5<-pchisq(surv_Pr_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_5<-pchisq(surv_Pr_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_5<-pchisq(surv_Pr_varied_primary_2.0_5$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_6<-pchisq(surv_Pr_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_6<-pchisq(surv_Pr_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_6<-pchisq(surv_Pr_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_6<-pchisq(surv_Pr_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_6<-pchisq(surv_Pr_varied_primary_2.0_6$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_7<-pchisq(surv_Pr_varied_primary_0.001_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_7<-pchisq(surv_Pr_varied_primary_0.01_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_7<-pchisq(surv_Pr_varied_primary_0.1_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_7<-pchisq(surv_Pr_varied_primary_1.0_7$chisq, df=1, lower.tail=F)
#pvalue_Pr_varied_primary_2.0_7<-pchisq(surv_Pr_varied_primary_2.0_7$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_8<-pchisq(surv_Pr_varied_primary_0.001_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_8<-pchisq(surv_Pr_varied_primary_0.01_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_8<-pchisq(surv_Pr_varied_primary_0.1_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_8<-pchisq(surv_Pr_varied_primary_1.0_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_8<-pchisq(surv_Pr_varied_primary_2.0_8$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.001_9<-pchisq(surv_Pr_varied_primary_0.001_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.01_9<-pchisq(surv_Pr_varied_primary_0.01_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_0.1_9<-pchisq(surv_Pr_varied_primary_0.1_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_1.0_9<-pchisq(surv_Pr_varied_primary_1.0_9$chisq, df=1, lower.tail=F)
pvalue_Pr_varied_primary_2.0_9<-pchisq(surv_Pr_varied_primary_2.0_9$chisq, df=1, lower.tail=F)# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Pr_varied_primary_0.001 <- c(pvalue_Pr_varied_primary_0.001_2 ,pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
test.stat_Pr_varied_primary_0.001 <- -2*sum(log(pvalues_Pr_varied_primary_0.001))
pfinal_Pr_varied_primary_0.001<-pchisq(test.stat_Pr_varied_primary_0.001, df=2*length(pvalues_Pr_varied_primary_0.001), lower.tail=F)
pvalues_Pr_varied_primary_0.01 <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
test.stat_Pr_varied_primary_0.01 <- -2*sum(log(pvalues_Pr_varied_primary_0.01))
pfinal_Pr_varied_primary_0.01<-pchisq(test.stat_Pr_varied_primary_0.01, df=2*length(pvalues_Pr_varied_primary_0.01), lower.tail=F)
pvalues_Pr_varied_primary_0.1 <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2,pvalue_Pr_varied_primary_0.1_3,pvalue_Pr_varied_primary_0.1_4,pvalue_Pr_varied_primary_0.1_5,pvalue_Pr_varied_primary_0.1_6,pvalue_Pr_varied_primary_0.1_7,pvalue_Pr_varied_primary_0.1_8,pvalue_Pr_varied_primary_0.1_9)
test.stat_Pr_varied_primary_0.1 <- -2*sum(log(pvalues_Pr_varied_primary_0.1))
pfinal_Pr_varied_primary_0.1<-pchisq(test.stat_Pr_varied_primary_0.1, df=2*length(pvalues_Pr_varied_primary_0.1), lower.tail=F)
pvalues_Pr_varied_primary_1.0 <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2,pvalue_Pr_varied_primary_1.0_3,pvalue_Pr_varied_primary_1.0_4,pvalue_Pr_varied_primary_1.0_5,pvalue_Pr_varied_primary_1.0_6,pvalue_Pr_varied_primary_1.0_7,pvalue_Pr_varied_primary_1.0_8,pvalue_Pr_varied_primary_1.0_9)
test.stat_Pr_varied_primary_1.0 <- -2*sum(log(pvalues_Pr_varied_primary_1.0))
pfinal_Pr_varied_primary_1.0<-pchisq(test.stat_Pr_varied_primary_1.0, df=2*length(pvalues_Pr_varied_primary_1.0), lower.tail=F)
pvalues_Pr_varied_primary_2.0 <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2,pvalue_Pr_varied_primary_2.0_3,pvalue_Pr_varied_primary_2.0_4,pvalue_Pr_varied_primary_2.0_5,pvalue_Pr_varied_primary_2.0_6,pvalue_Pr_varied_primary_2.0_8,pvalue_Pr_varied_primary_2.0_9)
test.stat_Pr_varied_primary_2.0 <- -2*sum(log(pvalues_Pr_varied_primary_2.0))
pfinal_Pr_varied_primary_2.0<-pchisq(test.stat_Pr_varied_primary_2.0, df=2*length(pvalues_Pr_varied_primary_2.0), lower.tail=F)pvalues_Pr_varied_primary_0.001_table <- c(pvalue_Pr_varied_primary_0.001_1, pvalue_Pr_varied_primary_0.001_2, pvalue_Pr_varied_primary_0.001_3, pvalue_Pr_varied_primary_0.001_4, pvalue_Pr_varied_primary_0.001_5, pvalue_Pr_varied_primary_0.001_6, pvalue_Pr_varied_primary_0.001_7, pvalue_Pr_varied_primary_0.001_8, pvalue_Pr_varied_primary_0.001_9)
p_Pr_varied_primary_0.001<-scales::pvalue(pvalues_Pr_varied_primary_0.001_table)
pvalues_Pr_varied_primary_0.01_table <- c(pvalue_Pr_varied_primary_0.01_1,pvalue_Pr_varied_primary_0.01_2,pvalue_Pr_varied_primary_0.01_3,pvalue_Pr_varied_primary_0.01_4,pvalue_Pr_varied_primary_0.01_5,pvalue_Pr_varied_primary_0.01_6,pvalue_Pr_varied_primary_0.01_7,pvalue_Pr_varied_primary_0.01_8,pvalue_Pr_varied_primary_0.01_9)
p_Pr_varied_primary_0.01<-scales::pvalue(pvalues_Pr_varied_primary_0.01_table)
pvalues_Pr_varied_primary_0.1_table <- c(pvalue_Pr_varied_primary_0.1_1,pvalue_Pr_varied_primary_0.1_2, pvalue_Pr_varied_primary_0.1_3, pvalue_Pr_varied_primary_0.1_4, pvalue_Pr_varied_primary_0.1_5, pvalue_Pr_varied_primary_0.1_6, pvalue_Pr_varied_primary_0.1_7, pvalue_Pr_varied_primary_0.1_8, pvalue_Pr_varied_primary_0.1_9)
p_Pr_varied_primary_0.1<-scales::pvalue(pvalues_Pr_varied_primary_0.1_table)
pvalues_Pr_varied_primary_1.0_table <- c(pvalue_Pr_varied_primary_1.0_1,pvalue_Pr_varied_primary_1.0_2, pvalue_Pr_varied_primary_1.0_3, pvalue_Pr_varied_primary_1.0_4, pvalue_Pr_varied_primary_1.0_5, pvalue_Pr_varied_primary_1.0_6, pvalue_Pr_varied_primary_1.0_7, pvalue_Pr_varied_primary_1.0_8, pvalue_Pr_varied_primary_1.0_9)
p_Pr_varied_primary_1.0<-scales::pvalue(pvalues_Pr_varied_primary_1.0_table)
pvalues_Pr_varied_primary_2.0_table <- c(pvalue_Pr_varied_primary_2.0_1,pvalue_Pr_varied_primary_2.0_2, pvalue_Pr_varied_primary_2.0_3, pvalue_Pr_varied_primary_2.0_4, pvalue_Pr_varied_primary_2.0_5, pvalue_Pr_varied_primary_2.0_6, NA, pvalue_Pr_varied_primary_2.0_8, pvalue_Pr_varied_primary_2.0_9)
p_Pr_varied_primary_2.0<-scales::pvalue(pvalues_Pr_varied_primary_2.0_table)
p_Pr_varied_primary_2.0 <- p_Pr_varied_primary_2.0 %>% replace_na("--")
pvalue_table_Pr_varied_primary <- data.frame("0.001" = p_Pr_varied_primary_0.001, "0.01" = p_Pr_varied_primary_0.01, "0.1" = p_Pr_varied_primary_0.1, "1.0" = p_Pr_varied_primary_1.0, "2.0" = p_Pr_varied_primary_2.0)
colnames(pvalue_table_Pr_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)
pfisher_Pr_varied_primary<- c(pfinal_Pr_varied_primary_0.001, pfinal_Pr_varied_primary_0.01, pfinal_Pr_varied_primary_0.1, pfinal_Pr_varied_primary_1.0, NA)
pfisher_Pr_varied_primary<-scales::pvalue(pfisher_Pr_varied_primary)
pfisher_Pr_varied_primary<- pfisher_Pr_varied_primary%>% replace_na("N.C.")
pvalue_table_Pr_varied_primary[nrow(pvalue_table_Pr_varied_primary) + 1,]<-pfisher_Pr_varied_primary
row.names(pvalue_table_Pr_varied_primary) <- c("Block 1 - October 18th, 2020", "Block 2 - October 20th, 2020", "Block 3 - October 29th, 2020","Block 4 - November 2nd, 2020", "Block 5 - November 7th, 2020", "Block 6 - November 28th, 2020","Block 7 - December 4th, 2020", "Block 8 - June 17th, 2021", "Block 9 - June 18th, 2021", " Final P-value ")
pvalue_table_Pr_varied_primary %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
pack_rows("Individual Blocks", 1, 9) %>%
pack_rows("Fisher's Combined", 10, 10)|
Primary Infectious Dose (Abs600nm)
|
|||||
|---|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | 2 | |
| Individual Blocks | |||||
| Block 1 - October 18th, 2020 | 0.072 | 0.293 | <0.001 | 0.070 | 0.041 |
| Block 2 - October 20th, 2020 | 0.526 | 0.211 | <0.001 | 0.002 | 0.072 |
| Block 3 - October 29th, 2020 | 0.628 | 0.018 | <0.001 | <0.001 | 0.012 |
| Block 4 - November 2nd, 2020 | 0.942 | 0.009 | <0.001 | <0.001 | 0.407 |
| Block 5 - November 7th, 2020 | 0.693 | 0.332 | <0.001 | <0.001 | 0.010 |
| Block 6 - November 28th, 2020 | 0.064 | 0.913 | <0.001 | <0.001 | 0.033 |
| Block 7 - December 4th, 2020 | 0.968 | 0.007 | <0.001 | <0.001 | – |
| Block 8 - June 17th, 2021 | 0.726 | 0.063 | <0.001 | <0.001 | <0.001 |
| Block 9 - June 18th, 2021 | 0.238 | 0.001 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | |||||
| Final P-value | 0.733 | <0.001 | <0.001 | <0.001 | N.C. |
- Conclusion: Moderate doses of the primary infection (Abs600nm = 0.01, 0.1 & 1.0) all provided significant protection against infection (P<0.001). The lowest dose (Abs600nm = 0.001) did not provide significant protection (p=0.76). The highest dose provided inconsistent protection with 8 out of the 9 replicates demonstrating better survival than controls and one with worse survival than controls therefore we could not calculate a combine p-value.
P. sneebia secondary infection
- Conclusion: All doses used to induce primary infection provided significant protection (P < 0.001), except for the lowest dose (Abs600nm = 0.261). Note: the highest primary dose provided inconsistent protection (4 out of 5 blocks had significant protection for all others, for Block 1 the mortality was the same as flies carrying only secondary infection)
Figure 5B
Ps_varied_primary <- subset(Ps_varied_primary, secondary == "P.snee")
table(Ps_varied_primary$date, Ps_varied_primary$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 3/10/2021 38 49 48 49 30 6
## 3/17/2021 37 50 46 42 22 0
## 3/31/2021 39 46 47 42 43 52
## 4/21/2021 35 44 46 45 30 49
## 4/28/2021 37 45 41 43 38 34
## 6/20/2021 27 48 49 48 47 60table(Ps_varied_primary$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 213 282 277 269 210 201#png(file="Figure5B.png", width = 1200, height = 900, res = 130)
par(mar=c(6,5,1,1)) #sets the margins around the graph so there is room for the labels
plot(survfit(Surv(Ps_varied_primary$hour, Ps_varied_primary$status) ~ Ps_varied_primary$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=5, yaxt='n', xaxt='n', ylab="proportion alive", xlab="hours post-secondary infection", cex.lab=2.5, xlim=c(0,48), xaxs='i', ylim=c(0, 1), yaxs='i')
axis(2, las=2, cex.axis=1.5, lwd=5)
axis(1, cex.axis=1.5, lwd=5)
box(lwd=5)
legend("bottomleft", bty="n" , c("secondary infection only", "chronic (0.001) & secondary infection", "chronic (0.01) & secondary infection", "chronic (0.1) & secondary infection", "chronic (1.0) & secondary infection", "chronic (2.0) & secondary infection"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1.5)
#dev.off()Individual Blocks
Block 1
Ps_varied_primary_block_1 <- subset(Ps_varied_primary, date == "3/10/2021")
table(Ps_varied_primary_block_1$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 38 49 48 49 30 6#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_1<-subset(Ps_varied_primary_block_1, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_1$primary_dose)##
## 0 0.001
## 38 49Ps_varied_primary_0.01_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_1<-subset(Ps_varied_primary_block_1, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_1<-subset(Ps_varied_primary_block_1, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_1)
surv_Ps_varied_primary_0.001_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 30.3 1.98 4.76
## primary=S.m 49 47 54.7 1.10 4.76
##
## Chisq= 4.8 on 1 degrees of freedom, p= 0.03surv_Ps_varied_primary_0.01_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_1)
surv_Ps_varied_primary_0.01_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 31.6 1.29 3.14
## primary=S.m 48 48 54.4 0.75 3.14
##
## Chisq= 3.1 on 1 degrees of freedom, p= 0.08surv_Ps_varied_primary_0.1_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_1)
surv_Ps_varied_primary_0.1_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 22.1 11.40 22.7
## primary=S.m 49 49 64.9 3.89 22.7
##
## Chisq= 22.7 on 1 degrees of freedom, p= 2e-06surv_Ps_varied_primary_1.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_1)
surv_Ps_varied_primary_1.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 36.7 0.0435 0.148
## primary=S.m 30 30 31.3 0.0511 0.148
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.7surv_Ps_varied_primary_2.0_1<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_1)
surv_Ps_varied_primary_2.0_1## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 38 38 37.12 0.0208 0.234
## primary=S.m 6 5 5.88 0.1314 0.234
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6plot(survfit(Surv(Ps_varied_primary_block_1$hour, Ps_varied_primary_block_1$status) ~ Ps_varied_primary_block_1$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 2
Ps_varied_primary_block_2 <- subset(Ps_varied_primary, date == "3/17/2021")
table(Ps_varied_primary_block_2$primary_dose)##
## 0 0.001 0.01 0.1 1
## 37 50 46 42 22#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_2<-subset(Ps_varied_primary_block_2, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_2$primary_dose)##
## 0 0.001
## 37 50Ps_varied_primary_0.01_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_2<-subset(Ps_varied_primary_block_2, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="1"|primary_dose=="0")
#Ps_varied_primary_2.0_2<-subset(Ps_varied_primary_block_2, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_2)
surv_Ps_varied_primary_0.001_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 38.5 0.0568 0.16
## primary=S.m 50 49 47.5 0.0460 0.16
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7surv_Ps_varied_primary_0.01_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_2)
surv_Ps_varied_primary_0.01_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 30 1.638 3.92
## primary=S.m 46 44 51 0.963 3.92
##
## Chisq= 3.9 on 1 degrees of freedom, p= 0.05surv_Ps_varied_primary_0.1_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_2)
surv_Ps_varied_primary_0.1_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 20.7 12.89 25.7
## primary=S.m 42 32 48.3 5.52 25.7
##
## Chisq= 25.7 on 1 degrees of freedom, p= 4e-07surv_Ps_varied_primary_1.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_2)
surv_Ps_varied_primary_1.0_2## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 27.8 3.04 9.36
## primary=S.m 22 14 23.2 3.64 9.36
##
## Chisq= 9.4 on 1 degrees of freedom, p= 0.002#surv_Ps_varied_primary_2.0_2<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_2)
#surv_Ps_varied_primary_2.0_2
plot(survfit(Surv(Ps_varied_primary_block_2$hour, Ps_varied_primary_block_2$status) ~ Ps_varied_primary_block_2$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 3
Ps_varied_primary_block_3 <- subset(Ps_varied_primary, date == "3/31/2021")
table(Ps_varied_primary_block_3$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 39 46 47 42 43 52#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_3<-subset(Ps_varied_primary_block_3, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_3$primary_dose)##
## 0 0.001
## 39 46Ps_varied_primary_0.01_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_3<-subset(Ps_varied_primary_block_3, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_3<-subset(Ps_varied_primary_block_3, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_3)
surv_Ps_varied_primary_0.001_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 23 0.177 0.383
## primary=S.m 46 29 27 0.151 0.383
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5surv_Ps_varied_primary_0.01_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_3)
surv_Ps_varied_primary_0.01_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 20.4 0.0196 0.0399
## primary=S.m 47 26 26.6 0.0150 0.0399
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8surv_Ps_varied_primary_0.1_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_3)
surv_Ps_varied_primary_0.1_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 23.3 0.222 0.504
## primary=S.m 42 30 27.7 0.186 0.504
##
## Chisq= 0.5 on 1 degrees of freedom, p= 0.5surv_Ps_varied_primary_1.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_3)
surv_Ps_varied_primary_1.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 12.3 6.08 11.3
## primary=S.m 43 8 16.7 4.50 11.3
##
## Chisq= 11.3 on 1 degrees of freedom, p= 8e-04surv_Ps_varied_primary_2.0_3<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_3)
surv_Ps_varied_primary_2.0_3## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 39 21 10.7 9.96 17.2
## primary=S.m 52 7 17.3 6.15 17.2
##
## Chisq= 17.2 on 1 degrees of freedom, p= 3e-05plot(survfit(Surv(Ps_varied_primary_block_3$hour, Ps_varied_primary_block_3$status) ~ Ps_varied_primary_block_3$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 4
Ps_varied_primary_block_4 <- subset(Ps_varied_primary, date == "4/21/2021")
table(Ps_varied_primary_block_4$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 35 44 46 45 30 49#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_4<-subset(Ps_varied_primary_block_4, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_4$primary_dose)##
## 0 0.001
## 35 44Ps_varied_primary_0.01_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_4<-subset(Ps_varied_primary_block_4, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_4<-subset(Ps_varied_primary_block_4, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_4)
surv_Ps_varied_primary_0.001_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 26.4 2.83 6.39
## primary=S.m 44 44 52.6 1.42 6.39
##
## Chisq= 6.4 on 1 degrees of freedom, p= 0.01surv_Ps_varied_primary_0.01_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_4)
surv_Ps_varied_primary_0.01_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 18.9 13.74 26
## primary=S.m 46 39 55.1 4.71 26
##
## Chisq= 26 on 1 degrees of freedom, p= 3e-07surv_Ps_varied_primary_0.1_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_4)
surv_Ps_varied_primary_0.1_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 12.4 41.3 69.8
## primary=S.m 45 22 44.6 11.5 69.8
##
## Chisq= 69.8 on 1 degrees of freedom, p= <2e-16surv_Ps_varied_primary_1.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_4)
surv_Ps_varied_primary_1.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 16 22.5 46.5
## primary=S.m 30 7 26 13.9 46.5
##
## Chisq= 46.5 on 1 degrees of freedom, p= 9e-12surv_Ps_varied_primary_2.0_4<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_4)
surv_Ps_varied_primary_2.0_4## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_4)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 35 35 19.5 12.33 23
## primary=S.m 49 26 41.5 5.79 23
##
## Chisq= 23 on 1 degrees of freedom, p= 2e-06plot(survfit(Surv(Ps_varied_primary_block_4$hour, Ps_varied_primary_block_4$status) ~ Ps_varied_primary_block_4$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 5
Ps_varied_primary_block_5 <- subset(Ps_varied_primary, date == "4/28/2021")
table(Ps_varied_primary_block_5$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 37 45 41 43 38 34#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_5<-subset(Ps_varied_primary_block_5, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_5$primary_dose)##
## 0 0.001
## 37 45Ps_varied_primary_0.01_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_5<-subset(Ps_varied_primary_block_5, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_5<-subset(Ps_varied_primary_block_5, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_5)
surv_Ps_varied_primary_0.001_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 37.2 0.000842 0.00276
## primary=S.m 45 45 44.8 0.000699 0.00276
##
## Chisq= 0 on 1 degrees of freedom, p= 1surv_Ps_varied_primary_0.01_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_5)
surv_Ps_varied_primary_0.01_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 36.4 0.00838 0.0276
## primary=S.m 41 41 41.6 0.00735 0.0276
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9surv_Ps_varied_primary_0.1_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_5)
surv_Ps_varied_primary_0.1_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 28.3 2.68 6.9
## primary=S.m 43 43 51.7 1.47 6.9
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.009surv_Ps_varied_primary_1.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_5)
surv_Ps_varied_primary_1.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 17.9 20.24 40.2
## primary=S.m 38 34 53.1 6.85 40.2
##
## Chisq= 40.2 on 1 degrees of freedom, p= 2e-10surv_Ps_varied_primary_2.0_5<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_5)
surv_Ps_varied_primary_2.0_5## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 37 37 18.8 17.71 37.4
## primary=S.m 34 30 48.2 6.89 37.4
##
## Chisq= 37.4 on 1 degrees of freedom, p= 9e-10plot(survfit(Surv(Ps_varied_primary_block_5$hour, Ps_varied_primary_block_5$status) ~ Ps_varied_primary_block_5$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Block 6
Ps_varied_primary_block_6 <- subset(Ps_varied_primary, date == "6/20/2021")
table(Ps_varied_primary_block_6$primary_dose)##
## 0 0.001 0.01 0.1 1 2
## 27 48 49 48 47 60#subset by primary dose for pairwise comparison dose
Ps_varied_primary_0.001_6<-subset(Ps_varied_primary_block_6, primary_dose == "0.001"|primary_dose=="0")
table(Ps_varied_primary_0.001_6$primary_dose)##
## 0 0.001
## 27 48Ps_varied_primary_0.01_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.01"|primary_dose=="0")
Ps_varied_primary_0.1_6<-subset(Ps_varied_primary_block_6, primary_dose=="0.1"|primary_dose=="0")
Ps_varied_primary_1.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="1"|primary_dose=="0")
Ps_varied_primary_2.0_6<-subset(Ps_varied_primary_block_6, primary_dose=="2"|primary_dose=="0")
surv_Ps_varied_primary_0.001_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.001_6)
surv_Ps_varied_primary_0.001_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.001_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 24.9 0.1730 0.436
## primary=S.m 48 48 50.1 0.0861 0.436
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5surv_Ps_varied_primary_0.01_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.01_6)
surv_Ps_varied_primary_0.01_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.01_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 20.2 2.262 4.59
## primary=S.m 49 46 52.8 0.868 4.59
##
## Chisq= 4.6 on 1 degrees of freedom, p= 0.03surv_Ps_varied_primary_0.1_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_0.1_6)
surv_Ps_varied_primary_0.1_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_0.1_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 10.5 25.69 40.5
## primary=S.m 48 31 47.5 5.71 40.5
##
## Chisq= 40.5 on 1 degrees of freedom, p= 2e-10surv_Ps_varied_primary_1.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_1.0_6)
surv_Ps_varied_primary_1.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_1.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 7.82 47.01 71.3
## primary=S.m 47 19 38.18 9.63 71.3
##
## Chisq= 71.3 on 1 degrees of freedom, p= <2e-16surv_Ps_varied_primary_2.0_6<-survdiff(Surv(day, status) ~ primary, data=Ps_varied_primary_2.0_6)
surv_Ps_varied_primary_2.0_6## Call:
## survdiff(formula = Surv(day, status) ~ primary, data = Ps_varied_primary_2.0_6)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## primary=PBS 27 27 14.1 11.88 19
## primary=S.m 60 37 49.9 3.35 19
##
## Chisq= 19 on 1 degrees of freedom, p= 1e-05plot(survfit(Surv(Ps_varied_primary_block_6$hour, Ps_varied_primary_block_6$status) ~ Ps_varied_primary_block_6$condition), col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"), lty=c(1, 1, 1, 1, 1, 1), lwd=3.4, yaxt='n', xaxt='n', ylab="Proportion Alive", xlab="Hours Post-Secondary Infection", cex.lab=1.3, xlim=c(10,48), xaxs='i', ylim=c(0, 1.05), yaxs='i')
axis(2, las=2, cex.axis=1, lwd=4.5)
axis(1, cex.axis=1, lwd=4.5)
box(lwd=4.5)
legend("bottomleft", bty="n" , c("PBS-P.snee 2.0", "S.m 0.001-P.snee 2.0", "S.m 0.01-P.snee 2.0", "S.m 0.1-P.snee 2.0", "S.m 1.0-P.snee 2.0", "S.m 2.0-P.snee 2.0"),col=c("gray", "#B5E8FC", "#75D6FD", "#09B6FC", "#055df5", "#003694"),lty=c(1, 1, 1, 1, 1, 1),lwd=5, cex=1)
Statistical Tests
- In order test whether there was protection provided by chronic infection for each secondary infectious dose, p-values from the log rank test for each date was combined using the Fisher’s test (experiments on all dates showed higher survival in conditions with chronic infection). Final p-value is reported in the text.
# Grab p-values from each individual logrank test (results shown above in each block)
pvalue_Ps_varied_primary_0.001_1<-pchisq(surv_Ps_varied_primary_0.001_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_1<-pchisq(surv_Ps_varied_primary_0.01_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_1<-pchisq(surv_Ps_varied_primary_0.1_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_1<-pchisq(surv_Ps_varied_primary_1.0_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_1<-pchisq(surv_Ps_varied_primary_2.0_1$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_2<-pchisq(surv_Ps_varied_primary_0.001_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_2<-pchisq(surv_Ps_varied_primary_0.01_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_2<-pchisq(surv_Ps_varied_primary_0.1_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_2<-pchisq(surv_Ps_varied_primary_1.0_2$chisq, df=1, lower.tail=F)
#pvalue_Ps_varied_primary_2.0_2<-pchisq(surv_Ps_varied_primary_2.0_2$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_3<-pchisq(surv_Ps_varied_primary_0.001_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_3<-pchisq(surv_Ps_varied_primary_0.01_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_3<-pchisq(surv_Ps_varied_primary_0.1_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_3<-pchisq(surv_Ps_varied_primary_1.0_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_3<-pchisq(surv_Ps_varied_primary_2.0_3$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_4<-pchisq(surv_Ps_varied_primary_0.001_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_4<-pchisq(surv_Ps_varied_primary_0.01_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_4<-pchisq(surv_Ps_varied_primary_0.1_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_4<-pchisq(surv_Ps_varied_primary_1.0_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_4<-pchisq(surv_Ps_varied_primary_2.0_4$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_5<-pchisq(surv_Ps_varied_primary_0.001_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_5<-pchisq(surv_Ps_varied_primary_0.01_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_5<-pchisq(surv_Ps_varied_primary_0.1_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_5<-pchisq(surv_Ps_varied_primary_1.0_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_5<-pchisq(surv_Ps_varied_primary_2.0_5$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.001_6<-pchisq(surv_Ps_varied_primary_0.001_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.01_6<-pchisq(surv_Ps_varied_primary_0.01_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_0.1_6<-pchisq(surv_Ps_varied_primary_0.1_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_1.0_6<-pchisq(surv_Ps_varied_primary_1.0_6$chisq, df=1, lower.tail=F)
pvalue_Ps_varied_primary_2.0_6<-pchisq(surv_Ps_varied_primary_2.0_6$chisq, df=1, lower.tail=F)# Use all p-values from a given condition to get Fisher's combined test p-value
pvalues_Ps_varied_primary_0.001 <- c(pvalue_Ps_varied_primary_0.001_2 ,pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
test.stat_Ps_varied_primary_0.001 <- -2*sum(log(pvalues_Ps_varied_primary_0.001))
pfinal_Ps_varied_primary_0.001<-pchisq(test.stat_Ps_varied_primary_0.001, df=2*length(pvalues_Ps_varied_primary_0.001), lower.tail=F)
pvalues_Ps_varied_primary_0.01 <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
test.stat_Ps_varied_primary_0.01 <- -2*sum(log(pvalues_Ps_varied_primary_0.01))
pfinal_Ps_varied_primary_0.01<-pchisq(test.stat_Ps_varied_primary_0.01, df=2*length(pvalues_Ps_varied_primary_0.01), lower.tail=F)
pvalues_Ps_varied_primary_0.1 <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2,pvalue_Ps_varied_primary_0.1_3,pvalue_Ps_varied_primary_0.1_4,pvalue_Ps_varied_primary_0.1_5,pvalue_Ps_varied_primary_0.1_6)
test.stat_Ps_varied_primary_0.1 <- -2*sum(log(pvalues_Ps_varied_primary_0.1))
pfinal_Ps_varied_primary_0.1<-pchisq(test.stat_Ps_varied_primary_0.1, df=2*length(pvalues_Ps_varied_primary_0.1), lower.tail=F)
pvalues_Ps_varied_primary_1.0 <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2,pvalue_Ps_varied_primary_1.0_3,pvalue_Ps_varied_primary_1.0_4,pvalue_Ps_varied_primary_1.0_5,pvalue_Ps_varied_primary_1.0_6)
test.stat_Ps_varied_primary_1.0 <- -2*sum(log(pvalues_Ps_varied_primary_1.0))
pfinal_Ps_varied_primary_1.0<-pchisq(test.stat_Ps_varied_primary_1.0, df=2*length(pvalues_Ps_varied_primary_1.0), lower.tail=F)
pvalues_Ps_varied_primary_2.0 <- c(pvalue_Ps_varied_primary_2.0_1,pvalue_Ps_varied_primary_2.0_3,pvalue_Ps_varied_primary_2.0_4,pvalue_Ps_varied_primary_2.0_5,pvalue_Ps_varied_primary_2.0_6)
test.stat_Ps_varied_primary_2.0 <- -2*sum(log(pvalues_Ps_varied_primary_2.0))
pfinal_Ps_varied_primary_2.0<-pchisq(test.stat_Ps_varied_primary_2.0, df=2*length(pvalues_Ps_varied_primary_2.0), lower.tail=F)pvalues_Ps_varied_primary_0.001_table <- c(pvalue_Ps_varied_primary_0.001_1, pvalue_Ps_varied_primary_0.001_2, pvalue_Ps_varied_primary_0.001_3, pvalue_Ps_varied_primary_0.001_4, pvalue_Ps_varied_primary_0.001_5, pvalue_Ps_varied_primary_0.001_6)
p_Ps_varied_primary_0.001<-scales::pvalue(pvalues_Ps_varied_primary_0.001_table)
pvalues_Ps_varied_primary_0.01_table <- c(pvalue_Ps_varied_primary_0.01_1,pvalue_Ps_varied_primary_0.01_2,pvalue_Ps_varied_primary_0.01_3,pvalue_Ps_varied_primary_0.01_4,pvalue_Ps_varied_primary_0.01_5,pvalue_Ps_varied_primary_0.01_6)
p_Ps_varied_primary_0.01<-scales::pvalue(pvalues_Ps_varied_primary_0.01_table)
pvalues_Ps_varied_primary_0.1_table <- c(pvalue_Ps_varied_primary_0.1_1,pvalue_Ps_varied_primary_0.1_2, pvalue_Ps_varied_primary_0.1_3, pvalue_Ps_varied_primary_0.1_4, pvalue_Ps_varied_primary_0.1_5, pvalue_Ps_varied_primary_0.1_6)
p_Ps_varied_primary_0.1<-scales::pvalue(pvalues_Ps_varied_primary_0.1_table)
pvalues_Ps_varied_primary_1.0_table <- c(pvalue_Ps_varied_primary_1.0_1,pvalue_Ps_varied_primary_1.0_2, pvalue_Ps_varied_primary_1.0_3, pvalue_Ps_varied_primary_1.0_4, pvalue_Ps_varied_primary_1.0_5, pvalue_Ps_varied_primary_1.0_6)
p_Ps_varied_primary_1.0<-scales::pvalue(pvalues_Ps_varied_primary_1.0_table)
pvalues_Ps_varied_primary_2.0_table <- c(pvalue_Ps_varied_primary_2.0_1,NA, pvalue_Ps_varied_primary_2.0_3, pvalue_Ps_varied_primary_2.0_4, pvalue_Ps_varied_primary_2.0_5, pvalue_Ps_varied_primary_2.0_6)
p_Ps_varied_primary_2.0<-scales::pvalue(pvalues_Ps_varied_primary_2.0_table)
p_Ps_varied_primary_2.0 <- p_Ps_varied_primary_2.0 %>% replace_na("--")
pvalue_table_Ps_varied_primary <- data.frame("0.001" = p_Ps_varied_primary_0.001, "0.01" = p_Ps_varied_primary_0.01, "0.1" = p_Ps_varied_primary_0.1, "1.0" = p_Ps_varied_primary_1.0, "2.0" = p_Ps_varied_primary_2.0)
colnames(pvalue_table_Ps_varied_primary) <- c(0.001,0.01, 0.1, 1.0,2.0)
pfisher_Ps_varied_primary<- c(pfinal_Ps_varied_primary_0.001, pfinal_Ps_varied_primary_0.01, pfinal_Ps_varied_primary_0.1, pfinal_Ps_varied_primary_1.0, pfinal_Ps_varied_primary_2.0)
pfisher_Ps_varied_primary<-scales::pvalue(pfisher_Ps_varied_primary)
pvalue_table_Ps_varied_primary[nrow(pvalue_table_Ps_varied_primary) + 1,]<-pfisher_Ps_varied_primary
row.names(pvalue_table_Ps_varied_primary) <- c("Block 1 - March 10th, 2021", "Block 2 - March 17th, 2021", "Block 3 - March 31st, 2021","Block 4 - April 21t, 2021", "Block 5 - April 28th, 2021", "Block 6 - June 20th, 2021", " Final P-value ")
pvalue_table_Ps_varied_primary %>%
kbl(caption = "<center><strong>P-values for log-rank tests on individual doses within each experiment<center><strong>") %>%
kable_classic(full_width = F, html_font = "Cambria", position = "center")%>%
add_header_above(c(" "=1, "Primary Infectious Dose (Abs600nm)" = 5)) %>%
pack_rows("Individual Blocks", 1, 6) %>%
pack_rows("Fisher's Combined", 7, 7)|
Primary Infectious Dose (Abs600nm)
|
|||||
|---|---|---|---|---|---|
| 0.001 | 0.01 | 0.1 | 1 | 2 | |
| Individual Blocks | |||||
| Block 1 - March 10th, 2021 | 0.029 | 0.076 | <0.001 | 0.700 | 0.628 |
| Block 2 - March 17th, 2021 | 0.689 | 0.048 | <0.001 | 0.002 | – |
| Block 3 - March 31st, 2021 | 0.536 | 0.842 | 0.478 | <0.001 | <0.001 |
| Block 4 - April 21t, 2021 | 0.011 | <0.001 | <0.001 | <0.001 | <0.001 |
| Block 5 - April 28th, 2021 | 0.958 | 0.868 | 0.009 | <0.001 | <0.001 |
| Block 6 - June 20th, 2021 | 0.509 | 0.032 | <0.001 | <0.001 | <0.001 |
| Fisher’s Combined | |||||
| Final P-value | 0.261 | <0.001 | <0.001 | <0.001 | <0.001 |
- Conclusion: All doses used to induce primary infection provided significant protection (P < 0.001), except for the lowest dose (Abs600nm = 0.01)
Figure 5
q <- ggdraw() +
draw_image("Figure5A.png")
r <- ggdraw() +
draw_image("Figure5B.png")
multiplot5 <- q/r
v<-multiplot5
#plot_annotation(tag_levels = 'A')
ggsave("Fig.5.Primary.pdf", v, device = "pdf")
v